41 CTAB on AU 111 21x12x2 LAMMPS

This notebook contains a set of code snippets and functions, supporting analysis of LAMMPS output. Exemplarily, the minimization, equilibration and production trajectories of one CTA+ ion in the vicinity of a 111 gold layer are investigated.

With the help of numpy, pandas and matplotlib, LAMMPS output energy conbtributions are evaluated and plotted. With the help of ase, parmed, nglview and ipywidgets, trajectories are visulized.

With the help of ase and asap, radial distribution functions, distances, displacements and diffusivities are evaluated.

In [1]:
## preferred installation method for netcdf on NEMO locally:
# module load mpi/openmpi/2.1-gnu-5.2
# module load {...}
# export CPPFLAGS="${CPPFLAGS} -I${MPI_INC_DIR}"
# export LDFLAGS="${LDFLAGS} -L${MPI_LIB_DIR}"
# pip install --user netCDF4

## alternatively:
# pip install --user --global-option=build_ext --global-option="-L${MPI_INC_DIR}" netCDF4
In [2]:
%lsmagic # notebook commands
Out[2]:
Available line magics:
%alias  %alias_magic  %autocall  %automagic  %autosave  %bookmark  %cat  %cd  %clear  %colors  %config  %connect_info  %cp  %debug  %dhist  %dirs  %doctest_mode  %ed  %edit  %env  %gui  %hist  %history  %killbgscripts  %ldir  %less  %lf  %lk  %ll  %load  %load_ext  %loadpy  %logoff  %logon  %logstart  %logstate  %logstop  %ls  %lsmagic  %lx  %macro  %magic  %man  %matplotlib  %mkdir  %more  %mv  %notebook  %page  %pastebin  %pdb  %pdef  %pdoc  %pfile  %pinfo  %pinfo2  %popd  %pprint  %precision  %profile  %prun  %psearch  %psource  %pushd  %pwd  %pycat  %pylab  %qtconsole  %quickref  %recall  %rehashx  %reload_ext  %rep  %rerun  %reset  %reset_selective  %rm  %rmdir  %run  %save  %sc  %set_env  %store  %sx  %system  %tb  %time  %timeit  %unalias  %unload_ext  %who  %who_ls  %whos  %xdel  %xmode

Available cell magics:
%%!  %%HTML  %%SVG  %%bash  %%capture  %%debug  %%file  %%html  %%javascript  %%js  %%latex  %%markdown  %%perl  %%prun  %%pypy  %%python  %%python2  %%python3  %%ruby  %%script  %%sh  %%svg  %%sx  %%system  %%time  %%timeit  %%writefile

Automagic is ON, % prefix IS NOT needed for line magics.
In [3]:
# list environment variables, like bash's 'printenv'
%env
Out[3]:
{'BASH_FUNC_module()': '() {  eval $(/opt/bwhpc/common/admin/modules/module-wrapper/modulecmd bash $*)\n}',
 'CC': 'gcc',
 'CFLAGS': '-O3 -fPIC',
 'CLICOLOR': '1',
 'CMAKE_BIN_DIR': '/opt/bwhpc/common/devel/cmake/3.9.0/bin',
 'CMAKE_DOC_DIR': '/opt/bwhpc/common/devel/cmake/3.9.0/doc',
 'CMAKE_HOME': '/opt/bwhpc/common/devel/cmake/3.9.0',
 'CMAKE_MAN_DIR': '/opt/bwhpc/common/devel/cmake/3.9.0/man',
 'CMAKE_VERSION': '3.9.0',
 'CPATH': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/include',
 'CPPFLAGS': '-I/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/include',
 'CXX': 'g++',
 'CXXFLAGS': '-O3 -fPIC',
 'DISPLAY': 'localhost:16.0',
 'ENV': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17',
 'F77': 'gfortran',
 'F90': 'gfortran',
 'FC': 'gfortran',
 'FCFLAGS': '-O3 -fPIC',
 'GIT_PAGER': 'cat',
 'GMXBIN': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/bin',
 'GMXDATA': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/share/gromacs',
 'GMXLDLIB': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/lib64',
 'GMXLIB': '/home/fr/fr_fr/fr_jh1130/gromacs/2018.1/share/gromacs/top',
 'GMXMAN': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/share/man',
 'GMX_WS': '/work/ws/nemo/fr_jh1130-201708-0/jobs/gmxlab',
 'GNU_BIN_DIR': '/opt/bwhpc/common/compiler/gnu/5.2.0/bin',
 'GNU_HOME': '/opt/bwhpc/common/compiler/gnu/5.2.0',
 'GNU_LIB_DIR': '/opt/bwhpc/common/compiler/gnu/5.2.0/lib64',
 'GNU_MAN_DIR': '/opt/bwhpc/common/compiler/gnu/5.2.0/man',
 'GNU_VERSION': '5.2.0',
 'GROMACS_BIN_DIR': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/bin',
 'GROMACS_HOME': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2',
 'GROMACS_LIB_DIR': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/lib64',
 'GROMACS_MAN_DIR': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/share/man',
 'GROMACS_VERSION': '2018.1',
 'GROUP_WS': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0',
 'HISTCONTROL': 'ignoredups',
 'HISTSIZE': '1000',
 'HOME': '/home/fr/fr_fr/fr_jh1130',
 'HOSTNAME': 'login2.nemo.privat',
 'HWLOC_HIDE_ERRORS': '1',
 'JPY_PARENT_PID': '150507',
 'LANG': 'en_US.UTF-8',
 'LC_ADDRESS': 'de_DE.UTF-8',
 'LC_IDENTIFICATION': 'de_DE.UTF-8',
 'LC_MEASUREMENT': 'de_DE.UTF-8',
 'LC_MONETARY': 'de_DE.UTF-8',
 'LC_NAME': 'de_DE.UTF-8',
 'LC_NUMERIC': 'de_DE.UTF-8',
 'LC_PAPER': 'de_DE.UTF-8',
 'LC_TELEPHONE': 'de_DE.UTF-8',
 'LC_TIME': 'de_DE.UTF-8',
 'LDFLAGS': '-L/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/lib',
 'LD_INCLUDE_PATH': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/include',
 'LD_LIBRARY_PATH': '/home/fr/fr_fr/fr_jh1130/.local/lib:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/lib:/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/lib:/opt/bwhpc/common/compiler/gnu/5.2.0/lib64:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/lib64:/opt/bwhpc/common/devel/python/3.6.3.20171114/lib',
 'LESSOPEN': '||/usr/bin/lesspipe.sh %s',
 'LIBRARY_PATH': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/lib',
 'LMPBIN': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/bin',
 'LOADEDMODULES': 'devel/python/3.6.3:compiler/gnu/5.2:mpi/openmpi/2.1-gnu-5.2:gromacs/2018.1-gnu-5.2:devel/cmake/3.9:lammps/16Mar18-gnu-5.2-openmpi-2.1:ase/3.15.0:devel/python/3.6-local',
 'LOGNAME': 'fr_jh1130',
 'LS_COLORS': 'rs=0:di=01;34:ln=01;36:mh=00:pi=40;33:so=01;35:do=01;35:bd=40;33;01:cd=40;33;01:or=40;31;01:mi=01;05;37;41:su=37;41:sg=30;43:ca=30;41:tw=30;42:ow=34;42:st=37;44:ex=01;32:*.tar=01;31:*.tgz=01;31:*.arc=01;31:*.arj=01;31:*.taz=01;31:*.lha=01;31:*.lz4=01;31:*.lzh=01;31:*.lzma=01;31:*.tlz=01;31:*.txz=01;31:*.tzo=01;31:*.t7z=01;31:*.zip=01;31:*.z=01;31:*.Z=01;31:*.dz=01;31:*.gz=01;31:*.lrz=01;31:*.lz=01;31:*.lzo=01;31:*.xz=01;31:*.bz2=01;31:*.bz=01;31:*.tbz=01;31:*.tbz2=01;31:*.tz=01;31:*.deb=01;31:*.rpm=01;31:*.jar=01;31:*.war=01;31:*.ear=01;31:*.sar=01;31:*.rar=01;31:*.alz=01;31:*.ace=01;31:*.zoo=01;31:*.cpio=01;31:*.7z=01;31:*.rz=01;31:*.cab=01;31:*.jpg=01;35:*.jpeg=01;35:*.gif=01;35:*.bmp=01;35:*.pbm=01;35:*.pgm=01;35:*.ppm=01;35:*.tga=01;35:*.xbm=01;35:*.xpm=01;35:*.tif=01;35:*.tiff=01;35:*.png=01;35:*.svg=01;35:*.svgz=01;35:*.mng=01;35:*.pcx=01;35:*.mov=01;35:*.mpg=01;35:*.mpeg=01;35:*.m2v=01;35:*.mkv=01;35:*.webm=01;35:*.ogm=01;35:*.mp4=01;35:*.m4v=01;35:*.mp4v=01;35:*.vob=01;35:*.qt=01;35:*.nuv=01;35:*.wmv=01;35:*.asf=01;35:*.rm=01;35:*.rmvb=01;35:*.flc=01;35:*.avi=01;35:*.fli=01;35:*.flv=01;35:*.gl=01;35:*.dl=01;35:*.xcf=01;35:*.xwd=01;35:*.yuv=01;35:*.cgm=01;35:*.emf=01;35:*.axv=01;35:*.anx=01;35:*.ogv=01;35:*.ogx=01;35:*.aac=01;36:*.au=01;36:*.flac=01;36:*.mid=01;36:*.midi=01;36:*.mka=01;36:*.mp3=01;36:*.mpc=01;36:*.ogg=01;36:*.ra=01;36:*.wav=01;36:*.axa=01;36:*.oga=01;36:*.spx=01;36:*.xspf=01;36:',
 'MAIL': '/var/spool/mail/fr_jh1130',
 'MANPATH': '/home/fr/fr_fr/fr_jh1130/.local/share/man:/opt/bwhpc/common/devel/cmake/3.9.0/man:/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/share/man:/opt/bwhpc/common/compiler/gnu/5.2.0/man:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/share/man:/opt/bwhpc/common/devel/python/3.6.3.20171114/share/man:/usr/share/man:/usr/local/share/man',
 'MATHEMATICA_HOME': '/opt/bwhpc/fr/math/mathematica/10.0.1',
 'MOABHOMEDIR': '/opt/moab',
 'MODULEPATH': '/home/fr/fr_fr/fr_jh1130/modulefiles:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/modulefiles:/opt/bwhpc/modulefiles/applications/info:/opt/bwhpc/modulefiles/applications/fr:/opt/bwhpc/modulefiles/applications/common:/opt/bwhpc/modulefiles/libraries/info:/opt/bwhpc/modulefiles/libraries/fr:/opt/bwhpc/modulefiles/libraries/common:/opt/bwhpc/modulefiles/development/info:/opt/bwhpc/modulefiles/development/fr:/opt/bwhpc/modulefiles/development/common:/opt/bwhpc/modulefiles/obsolete/info:/opt/bwhpc/modulefiles/obsolete/fr:/opt/bwhpc/modulefiles/obsolete/common',
 'MODULERCFILE': '/opt/bwhpc/common/etc/modulerc',
 'MODULESHOME': '/usr/share/Modules',
 'MPIDIR': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2',
 'MPI_BIN_DIR': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/bin',
 'MPI_EXA_DIR': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/examples',
 'MPI_HOME': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2',
 'MPI_INC_DIR': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/include',
 'MPI_LIB_DIR': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/lib',
 'MPI_MAN_DIR': '/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/share/man',
 'MPI_VERSION': '2.1.3-gnu-5.2',
 'MPLBACKEND': 'module://ipykernel.pylab.backend_inline',
 'PAGER': 'cat',
 'PATH': '/opt/bwhpc/common/devel/python/3.6.3.20171114/bin:/home/fr/fr_fr/fr_jh1130/.local/bin:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/ase/3.15.0/:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/ase/3.15.0/bin:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/bin:/opt/bwhpc/common/devel/cmake/3.9.0/bin:/opt/bwhpc/common/mpi/openmpi/2.1.3-gnu-5.2/bin:/opt/bwhpc/common/compiler/gnu/5.2.0/bin:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/bin:/opt/bwhpc/common/devel/python/3.6.3.20171114/bin:/home/fr/fr_fr/fr_jh1130/bin:/opt/VirtualGL/bin:/usr/lib64/qt-3.3/bin:/opt/moab/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/usr/local/bin',
 'PERL5LIB': '/opt/moab/lib/perl5',
 'PKG_CONFIG_PATH': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/gromacs/2018.1-gnu-5.2/lib64/pkgconfig',
 'PREFIX': '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1',
 'PWD': '/work/ws/nemo/fr_jh1130-201708-0',
 'PYTHONPATH': '/home/fr/fr_fr/fr_jh1130/.local/lib/python3.6:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/ase/3.15.0/lib/python3.6/site-packages:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_modules/lammps/16Mar18-gnu-5.2-openmpi-2.1/lib/python3.6/site-packages:/opt/bwhpc/common/devel/python/3.6.3.20171114/lib/python3.6',
 'QTDIR': '/usr/lib64/qt-3.3',
 'QTINC': '/usr/lib64/qt-3.3/include',
 'QTLIB': '/usr/lib64/qt-3.3/lib',
 'QT_GRAPHICSSYSTEM_CHECKED': '1',
 'SHELL': '/bin/bash',
 'SHLVL': '2',
 'SSH_CLIENT': '2.202.161.150 57392 22',
 'SSH_CONNECTION': '2.202.161.150 57392 132.230.222.13 22',
 'SSH_TTY': '/dev/pts/5',
 'TERM': 'xterm-color',
 'USER': 'fr_jh1130',
 'WS': '/work/ws/nemo/fr_jh1130-201708-0',
 'XDG_RUNTIME_DIR': '/run/user/901576',
 'XDG_SESSION_ID': '203152',
 '_': '/opt/bwhpc/common/devel/python/3.6.3.20171114/bin/jupyter',
 '_LMFILES_': '/opt/bwhpc/modulefiles/development/common/devel/python/3.6.3:/opt/bwhpc/modulefiles/development/common/compiler/gnu/5.2:/opt/bwhpc/modulefiles/libraries/common/mpi/openmpi/2.1-gnu-5.2:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/modulefiles/gromacs/2018.1-gnu-5.2:/opt/bwhpc/modulefiles/development/common/devel/cmake/3.9:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/modulefiles/lammps/16Mar18-gnu-5.2-openmpi-2.1:/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/modulefiles/ase/3.15.0:/home/fr/fr_fr/fr_jh1130/modulefiles/devel/python/3.6-local',
 'python_BIN_DIR': '/opt/bwhpc/common/devel/python/3.6.3.20171114/bin',
 'python_HOME': '/opt/bwhpc/common/devel/python/3.6.3.20171114',
 'python_LDLIB_DIR': '/opt/bwhpc/common/devel/python/3.6.3.20171114/lib',
 'python_LIB_DIR': '/opt/bwhpc/common/devel/python/3.6.3.20171114/lib',
 'python_MAN_DIR': '/opt/bwhpc/common/devel/python/3.6.3.20171114/share/man',
 'python_VERSION': '3.6.3'}
In [4]:
%load_ext memory_profiler
In [5]:
# for some reason, nglview sometimes changes into some temporary directory
# therefore ALWAY use absolute filenames and paths
# %cd /work/ws/nemo/fr_jh1130-201708-0/jobs/lmplab/sds/201806/1_SDS_on_AU_100_1x4x4/

Imports

In [5]:
# system basics
import os
absolute_prefix = os.getcwd() # might be handy to get back to the initial working directory at any point
os.sep # '\' ond windows and '/' on unix-like
Out[5]:
'/'
In [6]:
# data analysis 
import pandas as pd 
import numpy as np

import ase
from asap3.analysis.rdf import RadialDistributionFunction

# file formats, input - output
import ase.io
from ase.io import read
from ase.io import NetCDFTrajectory
import parmed as pmd

# visualization
from ase.visualize import view
import nglview as nv
import matplotlib.pyplot as plt
import ipywidgets # just for jupyter notebooks
In [7]:
import ipyparallel as ipp

Global options

In [8]:
%matplotlib inline
In [9]:
# matplotlib settings

# expecially for presentation, larger font settings for plotting are recommendable
SMALL_SIZE = 14
MEDIUM_SIZE = 16
BIGGER_SIZE = 18

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure titlex

plt.rcParams["figure.figsize"] = (11,7) # the standard figure size
In [10]:
# numpy truncates the output of large array above the treshold length
np.set_printoptions(threshold=100) 

Definition of helper functions

In [11]:
# conversion units, only for better readability
fs = 1e-15 # s
ps = 1e-12 # s
In [12]:
AA = 1e-10 # m
In [13]:
def fullprint(*args, **kwargs):
    """prints a long numpy array without altering numpy treshold options permanently"""
    from pprint import pprint
    opt = np.get_printoptions()
    np.set_printoptions(edgeitems=3,infstr='inf', 
        linewidth=75, nanstr='nan', precision=8,
        suppress=False, threshold=100000, formatter=None)
    pprint(*args, **kwargs)
    np.set_printoptions(**opt)
In [14]:
def runningMeanFast(x, N):
    """a quick way to compute the running or rolling mean on a numpy array"""
    return np.convolve(x, np.ones((N,))/N)[(N-1):]
In [15]:
def running_mean(x, N):
    """another quick way to compute the running or rolling mean on a numpy array"""
    cumsum = np.cumsum(np.insert(x, 0, 0)) 
    return (cumsum[N:] - cumsum[:-N]) / float(N)
In [16]:
def subplotPosition(rows,cols):
    """generator for subplot positions"""
    for p in range(0,rows*cols):
        yield rows*100+cols*10+p+1

def addSubplot(x, y,
               title=None, xlabel=None, ylabel=None, legend=None,
               fig=None, ax=None, pos=None, figsize=(8,5)):
    """facilitate matplotlib figure & subplot creation. only one data series per call."""
    
    if not pos:
        pos = 111
    if not fig and not ax:
        fig = plt.figure(figsize=figsize)
    elif not fig and ax:
        fig = ax.get_figure()
    if not ax:
        ax = fig.add_subplot(pos)
      
    if legend:
        ax.plot(x,y,label=legend)
    else:
        ax.plot(x,y)
               
    if title:
        ax.set_title(title)
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
        
    if legend:
        ax.legend()
        
    return fig, ax
In [17]:
def makeThermoPlotsFromDataFrame(df, fig=None,
        time_label          = r'$\frac{\mathrm{Steps}}{2 \mathrm{fs}}$',
        temperature_label   = r'$\frac{T}{\mathrm{K}}$',
        pressure_label      = r'$\frac{P}{\mathrm{atm}}$',
        energy_label        = r'$\frac{E}{\mathrm{Kcal} \cdot \mathrm{mole}^{-1}}$'):
    """Automizes the plotting of thermo output."""

    rows = 3
    cols = 2
    if fig == None:
        fig = plt.figure(figsize=(cols*8,rows*5))
    
    def subplotPosition(rows,cols):
        for p in range(0,rows*cols):
            yield rows*100+cols*10+p+1
     
    def addSubplot(df,fig,pos,title,xlabel,ylabel):
        ax = fig.add_subplot(pos)
        df.plot(ax=ax) # taimed temperature
        ax.set_title(title)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
    
    pos = subplotPosition(rows,cols)
    
    # sum up intramolecular contributions 
    df["E_intramolecular"] = df[["E_bond","E_angle","E_dihed"]].sum(axis=1)
    
    addSubplot(df[["Temp"]],
               fig, next(pos), "Temperature", time_label, temperature_label)
    addSubplot(df[["Press"]],
               fig, next(pos), "Pressure", time_label, pressure_label)
   
    # intramolecular contributions (without angle)
    addSubplot(df[["E_intramolecular","E_bond","E_angle","E_dihed"]],
               fig, next(pos), "Intramolecular energies", time_label, energy_label)
    # intermolecular ("non-bonded") energy contribtutions
    # E_pair is the sum of the three latter, just as E_intramolecular in the plot above
    addSubplot(df[["E_pair","E_vdwl","E_coul","E_long"]],
                fig, next(pos), "Intermolecular (non-bonded) energies", time_label, energy_label)
    
    # visualize the difference between total and non-bonded potential:
    addSubplot(df[["PotEng","E_pair"]],
               fig, next(pos), "Total potential and non-bonded potential", time_label, energy_label)
    
    addSubplot(df[["TotEng","KinEng","PotEng"]],
               fig, next(pos), "Total, kinetic and potential energies", time_label, energy_label)
    
    fig.tight_layout()
    return fig
In [18]:
def makeRollingAverageThermoPlotsFromDataFrame(df, fig=None,
        time_label          = r'$\frac{\mathrm{Steps}}{2 \mathrm{fs}}$',
        temperature_label   = r'$\frac{T}{\mathrm{K}}$',
        pressure_label      = r'$\frac{P}{\mathrm{atm}}$',
        energy_label        = r'$\frac{E}{\mathrm{Kcal} \cdot \mathrm{mole}^{-1}}$',
        window = 1000):
    """Automizes the plotting of thermo output. Displays the rolling average with default window = 1000."""
    
    rows = 3
    cols = 2
    if fig == None:
        fig = plt.figure(figsize=(cols*8,rows*5))
    
    def subplotPosition(rows,cols):
        for p in range(0,rows*cols):
            yield rows*100+cols*10+p+1
     
    def addSubplot(df,fig,pos,title,xlabel,ylabel):
        ax = fig.add_subplot(pos)
        df.plot(ax=ax)
        ax.set_title(title)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
    
    pos = subplotPosition(rows,cols)
    
    # sum up intramolecular contributions 
    df["E_intramolecular"] = df[["E_bond","E_angle","E_dihed"]].sum(axis=1)
    
    addSubplot(df[["Temp"]].rolling(window=window,center=True).mean(),
               fig, next(pos), "Temperature", time_label, temperature_label)
    addSubplot(df[["Press"]].rolling(window=window,center=True).mean(),
               fig, next(pos), "Pressure", time_label, pressure_label)
   
    # intramolecular contributions (without angle)
    addSubplot(
      df[["E_intramolecular","E_bond","E_angle","E_dihed"]].rolling(window=window,center=True).mean(),
      fig, next(pos), "Intramolecular energies", time_label, energy_label)
    # intermolecular ("non-bonded") energy contribtutions
    # E_pair is the sum of the three latter, just as E_intramolecular in the plot above
    addSubplot(df[["E_pair","E_vdwl","E_coul","E_long"]].rolling(window=window,center=True).mean(),
                fig, next(pos), "Intermolecular (non-bonded) energies", time_label, energy_label)
    
    # visualize the difference between total and non-bonded potential:
    addSubplot(df[["PotEng","E_pair"]].rolling(window=window,center=True).mean(),
               fig, next(pos), "Total potential and non-bonded potential", time_label, energy_label)
    
    addSubplot(df[["TotEng","KinEng","PotEng"]].rolling(window=window,center=True).mean(),
               fig, next(pos), "Total, kinetic and potential energies", time_label, energy_label)
    
    fig.tight_layout() # tigh_layout avoids label overlap
    return fig
In [19]:
# ASE by default infers elements from LAMMPS atom types, in our case they are unrelated
# During preprocessing, our system went through several formats, one of them the 
# archaic .pdb format. Although ASE offers a pdb reader, it fails on our system.

# On the other hand, ParmEd is able to read .pdb and infer elements more or less accurately,
# but cannot process netCDF. Thus we combine both:
def inferTypes2NumbersFromPdbAndLmp(pdb_file,lmp_data_file):
    """Uses parmed's ability to infer elements from pdb files and constructs a dictionary and
    dictionary-like array for atom type -> atomic number assignments"""
    struct_pdb_pmd = pmd.read_PDB(pdb_file)
    struct_lmp_data = ase.io.read(lmp_data_file,format='lammps-data')
    
    resnames = np.unique([r.name for r in struct_pdb_pmd.residues])
    print("PDB contains following residue types {}".format(resnames))
    
    ions = [ a for a in struct_pdb_pmd.atoms if a.residue.name == 'ION' ] # sodium counterions
    
    # sodium correction: apparently, SOD and S are both interpreted as sulfur
    # bromide correction: in the pdb, the bromide ion is marked as CLA (chloride ion),
    #   because charmm2lammps.pl only supports placing Na+ or Cl- as counterions when preparing the system
    #   The chloride paramters have been changed to Horinek's bromide parameters manually after conversion 
    #   within the LAMMPS data file. What is more, the element C (carbon) is wrongly inferred 
    #   from residue name CLA by Parmed. Thus, atomic number 6 (C) is changed to Bromide her:
    for ion in ions:
        if ion.atomic_number == 16: # wrong: inferred S (sulfur)
            ion.atomic_number = 11 # Na
        elif ion.atomic_number == 6: # wrong: inferred C (sulfur)
            ion.atomic_number = 35 # Br
            
    # elements numbered as in periodic table
    atomic_numbers = np.array([ a.atomic_number for a in struct_pdb_pmd.atoms ]) 
    atomic_types = struct_lmp_data.get_atomic_numbers() # types as numbered in LAMMPS
    types2numbers = dict(zip(atomic_types,atomic_numbers)) # automatically picks unique tuples
    
    print("System contains {:d} atom types.".format(len(types2numbers)))
    
    # construct array, where indices represent LAMMPS type numbers and point to atomic numbers
    types2numbers_array = np.zeros(atomic_types.max()+1,dtype=np.uint)

    for k,v in types2numbers.items():
        types2numbers_array[k] = v
        
    # this kind of array representation allows for simple type2number conversion via an expression like
    #   types2numbers_array[atomic_types]
    
    return types2numbers, types2numbers_array
In [20]:
# helper function needs:
# trajectory, segement length, indices, element tuple list
# optional: start, end, rMax, nBins

def piecewiseRDF(traj, atom_indices, element_tuples,
                nSegment = 1000, nStart = 0, nEnd = None,
                rMax = 20, nBins = 1000):
    """Computes time-segment-wise element-element RDFs within a certain group of atoms of a trajectory"""

    # rMax is the rdf length (in Angstrom for LAMMPS output in real units)
    # nBins: can be understood as the number of data points on the RDF
    if not nEnd:
        nEnd = len(traj)-1

    nRDFs = len(element_tuples)

    listOfRdfLists = []
    for n in range(0,nRDFs):
        listOfRdfLists.append([])

    # actual distances
    rdf_x = (np.arange(nBins) + 0.5) * rMax / nBins

    # instead of computing an average rdf over the whole trajectory, 
    # we split the trajectory into several timespans of nSegement timestep length
    for curStart in range(nStart,nEnd,nSegment):
        print(curStart) # some progress report
        rdfObj = None
        for frame in traj[curStart:(curStart+nSegment)]:
            # the asap rdf functionality is not that convenient, but explicitely choosing
            # only the atoms we are interested in a priori, we can get exactly the rdf we want by 
            # making use of the "elements" option
            if rdfObj is None:
                rdfObj = RadialDistributionFunction(frame[atom_indices],
                                         rMax = rMax, nBins = nBins)
            else:
                rdfObj.atoms = frame[atom_indices]
            rdfObj.update()           

        # np.where facilitates the selection of according atom numbers by specifying the chemical symbol
        for n in range(0,nRDFs):
            curRdf = rdfObj.get_rdf(elements=element_tuples[n])
            listOfRdfLists[n].append(curRdf)

    return listOfRdfLists, rdf_x, rdfObj
In [21]:
# code snippet for neat plotting of all time-segemtn rdfs
def plotPiecewiceRdf(rdf_x, listOfRdfLists, legend=None,
                    nSegment = 1000, nStart=0, nEnd = None, cols = 2):
    """Plots a set of rdf constructed by piecewiseRDF(...)"""
    # Todo: Implement arbitrary start and end points

    
    N = 0
    for rdfList in listOfRdfLists:
        if len(rdfList) > N:
            N = len(rdfList)
            
    rows = np.ceil(N/cols).astype(int)
    
    pos = subplotPosition(rows,cols)
    fig = plt.figure(figsize=(rows*5,cols*8))

    for i in range(0,N):
        p = next(pos)
        ax = None
        for j, rdfList in enumerate(listOfRdfLists):
            if i < len(rdfList):
                curLegend = None
                if legend and j < len(legend):
                    curLegend = legend[j]
                    
                if not ax:
                    _, ax = addSubplot(rdf_x,rdfList[i], 
                               legend = curLegend, 
                               xlabel = r'$\frac{r}{\AA}$', 
                               ylabel='arbitrary density', 
                               title = "{} ps - {} ps".format(i*nSegment,(i+1)*nSegment),
                               fig = fig, pos = p)
                else:
                    _, _ = addSubplot(rdf_x, rdfList[i], 
                              legend = curLegend,
                              ax = ax, pos = p)

    fig.tight_layout()
    return fig
In [22]:
def piecewiseAveragedDistance(traj, reference_index, atom_indices,
                nSegment = 1000, nStart = 0, nEnd = None):
    """Piecewise averaged distance from one reference atom to a group of other atoms"""
    if not nEnd:
        nEnd = len(traj)-1

   
    N = (np.floor(nEnd - nStart)/nSegment).astype(int)
    
    t = (np.arange(nStart,nEnd,nSegment) + nSegment / 2)

    averageDistance = np.atleast_2d(np.zeros((3,N)))
         
    for i, curStart in enumerate(range(nStart,nEnd,nSegment)):
        #print(curStart) # some progress report
        curAverageDistance = np.atleast_2d(np.zeros((3,nSegment)))
        for j,f in enumerate(traj[curStart:(curStart+nSegment)]):
            curAverageDistance[:,j] = np.abs( f.get_distances(
                reference_index, atom_indices, mic=True, vector=True).mean(axis=0) )
        averageDistance[:,i] = curAverageDistance.mean(axis=1)
         
    return averageDistance, t
In [23]:
def piecewiseAveragedComComDistance(traj, group1_indices, group2_indices,
                nSegment = 1000, nStart = 0, nEnd = None):
    """Piecewise averaged distance from one atom group's center of mass to another one's"""

    if not nEnd:
        nEnd = len(traj)-1

   
    N = (np.floor(nEnd - nStart)/nSegment).astype(int)
    
    t = (np.arange(nStart,nEnd,nSegment) + nSegment / 2)

    averageDistance = np.atleast_2d(np.zeros((3,N)))
         
    for i, curStart in enumerate(range(nStart,nEnd,nSegment)):
        #print(curStart) # some progress report
        curAverageDistance = np.atleast_2d(np.zeros((3,nSegment)))
        for j,f in enumerate(traj[curStart:(curStart+nSegment)]):
            curAverageDistance[:,j] = np.abs(
                f[group1_indices].get_center_of_mass() - f[group2_indices].get_center_of_mass() )
        averageDistance[:,i] = curAverageDistance.mean(axis=1)
         
    return averageDistance, t
In [24]:
def comDisplacement(traj, atom_indices, dt = 10):
    """Evaluates the displacement of an atom group's COM 
    between each frame of the trajectory and another frame 'dt' indices later."""
    N = len(traj) - dt

    displacement = np.atleast_2d(np.zeros((3,N)))

    for j in range(N):
        reference_com   = traj[j][atom_indices].get_center_of_mass()
        dt_com           = traj[j+dt][atom_indices].get_center_of_mass()
        displacement[:,j] = np.abs(dt_com - reference_com)         
    
    return displacement
In [25]:
def evaluateDisplacement(displacement, dt = 10, window = 500):
    """Evaluates an anisotropic displacement vector and also computes diffusivities
    based on Enstein relation"""
    isotropic_displacement = np.linalg.norm(displacement,axis=0)

    time = np.arange( len(isotropic_displacement) )* TimePerFrame/ps

    rows = 3
    cols = 1
    
    fig = plt.figure(figsize=(8*cols, 5*rows))
    
    pos = subplotPosition(cols=cols,rows=rows)
    
    p = next(pos)
    _, ax = addSubplot(time[:-window+1], running_mean( isotropic_displacement,window), 
            title = "{:.2f} ps displacement".format(dt*TimePerFrame/ps), 
            xlabel = "time t / ps", ylabel= r'displacement $\frac{r}{\AA}$',
            legend = 'isotropic', fig = fig, pos = p)

    for i in range(0,3):
        addSubplot(time[:-window+1], 
                   running_mean(displacement[i,:], window), 
                   legend = distanceLabels[i], ax = ax, pos = p)
        
    p = next(pos)
    _, ax = addSubplot(time[:-window+1], running_mean( isotropic_displacement**2 / 3.0, window), 
            title = "{:.2f} ps MSD".format(dt*TimePerFrame/ps), 
            xlabel = "time t / ps", ylabel= r'MSD $\frac{r^2}{\AA^2}$',
            legend = 'isotropic', fig = fig, pos = p)
    
    for i in range(0,3):
        addSubplot(time[:-window+1], 
                   running_mean(displacement[i,:]**2, window), 
                   legend = distanceLabels[i], ax = ax, pos = p)
        
        
    # Einstein relation and converion to SI units
    D    = 1.0/2.0 * displacement**2 * AA**2 / (dt*TimePerFrame) 
    Diso = 1.0/6.0 * isotropic_displacement**2 * AA**2 / (dt*TimePerFrame) 
    
    p = next(pos)       
    _, ax = addSubplot(time[:-window+1], running_mean( Diso, window), 
            title = "diffusivities from {:.2f} ps MSD".format(dt*TimePerFrame/ps), 
            xlabel = "time t / ps", ylabel= r'D $\frac{m^2}{s}$',
            legend = 'isotropic', fig = fig, pos = p)

    
    for i in range(0,3):
        addSubplot(time[:-window+1], 
                   running_mean(D[i,:], window), 
                   legend = distanceLabels[i], ax = ax, pos = p)
        
    fig.tight_layout()
    
    return fig

Energy evaluations with pandas

In [58]:
dt = 2e-15 # s, 2 fs timestep

Minimization

In [59]:
%ls -lh *minimization.log
-rw-r--r-- 1 fr_jh1130 fr_fr 20K  7. Aug 16:10 41_CTAB_on_AU_111_21x12x2_bilayer_01_minimization.log
In [60]:
!./extract_thermo.sh *minimization.log 01_minimization_thermo.out
In [61]:
minimization_thermo_file = absolute_prefix + os.sep + '01_minimization_thermo.out'
In [62]:
minimization_thermo_pd = pd.read_csv(minimization_thermo_file,delim_whitespace=True)
In [63]:
minimization_thermo_pd.set_index("Step",inplace=True)
In [64]:
minimization_thermo_pd
Out[64]:
TotEng KinEng PotEng Temp Press Enthalpy E_bond E_angle E_dihed E_impro E_pair E_vdwl E_coul E_long E_tail Volume
Step
0 -195771.84 0 -195771.84 0 -1367.2245 -208799.36 31369.7640 8038.1440 337.88654 0 -235517.63 -27510.412 799137.94 -1007145.2 0 653350.57
100 -344887.99 0 -344887.99 0 -6936.7249 -410984.21 7829.3304 4814.3542 376.11682 0 -357907.79 -94068.576 753174.37 -1017013.6 0 653350.57
200 -356271.21 0 -356271.21 0 -5552.2242 -409175.29 8497.1175 5438.4882 375.12159 0 -370581.94 -87780.494 734973.10 -1017774.5 0 653350.57
300 -359738.20 0 -359738.20 0 -5241.3498 -409680.13 8708.1290 5559.3422 372.02656 0 -374377.70 -86479.091 730048.46 -1017947.1 0 653350.57
400 -361306.95 0 -361306.95 0 -5141.8331 -410300.64 8825.7396 5595.6140 369.51019 0 -376097.81 -85779.020 727684.86 -1018003.6 0 653350.57
500 -362138.47 0 -362138.47 0 -5097.1431 -410706.33 8870.3478 5631.2802 370.83873 0 -377010.94 -85548.472 726575.25 -1018037.7 0 653350.57
600 -362661.86 0 -362661.86 0 -5062.2305 -410897.06 8900.6397 5649.4127 371.60024 0 -377583.51 -85388.995 725859.67 -1018054.2 0 653350.57
700 -362992.87 0 -362992.87 0 -5049.7683 -411109.33 8915.4882 5647.2942 371.63893 0 -377927.29 -85295.070 725432.74 -1018065.0 0 653350.57
800 -363224.93 0 -363224.93 0 -5004.1837 -410907.04 8938.8342 5664.3778 372.98790 0 -378201.13 -85123.503 724995.63 -1018073.3 0 653350.57
900 -363392.10 0 -363392.10 0 -4989.4232 -410933.56 8943.8120 5669.4326 373.28253 0 -378378.63 -85069.695 724773.12 -1018082.1 0 653350.57
1000 -363478.02 0 -363478.02 0 -4980.3869 -410933.38 8951.4236 5670.8826 373.30451 0 -378473.64 -85036.204 724647.52 -1018084.9 0 653350.57
1100 -363538.79 0 -363538.79 0 -4964.1476 -410839.42 8956.0291 5672.6905 374.41680 0 -378541.93 -84988.725 724530.75 -1018084.0 0 653350.57
1200 -363588.08 0 -363588.08 0 -4977.4401 -411015.36 8955.9108 5673.0495 374.80939 0 -378591.85 -85026.607 724519.77 -1018085.0 0 653350.57
1300 -363636.53 0 -363636.53 0 -4974.3947 -411034.79 8958.4669 5676.1403 375.16095 0 -378646.30 -85016.954 724457.57 -1018086.9 0 653350.57
1400 -363679.11 0 -363679.11 0 -4957.7946 -410919.20 8962.4910 5676.6726 374.54767 0 -378692.82 -84976.638 724372.22 -1018088.4 0 653350.57
1500 -363698.32 0 -363698.32 0 -4950.9290 -410872.99 8964.0661 5677.7828 373.83962 0 -378714.01 -84958.783 724334.02 -1018089.3 0 653350.57
1600 -363721.37 0 -363721.37 0 -4935.3870 -410747.95 8966.9880 5680.2275 373.53519 0 -378742.12 -84920.331 724267.03 -1018088.8 0 653350.57
1700 -363740.86 0 -363740.86 0 -4933.5413 -410749.85 8967.6124 5680.6921 373.93221 0 -378763.10 -84915.910 724242.54 -1018089.7 0 653350.57
1800 -363759.75 0 -363759.75 0 -4918.8334 -410628.60 8970.5268 5683.2477 374.17921 0 -378787.71 -84875.284 724177.92 -1018090.3 0 653350.57
1900 -363780.76 0 -363780.76 0 -4905.3757 -410521.37 8973.8424 5684.1327 374.18647 0 -378812.92 -84835.674 724113.41 -1018090.6 0 653350.57
2000 -363805.66 0 -363805.66 0 -4907.4203 -410565.76 8973.9439 5686.4010 375.04510 0 -378841.05 -84846.750 724096.74 -1018091.0 0 653350.57
2100 -363823.81 0 -363823.81 0 -4908.0472 -410589.88 8975.3481 5687.5704 375.27501 0 -378862.00 -84837.350 724066.16 -1018090.8 0 653350.57
2200 -363835.80 0 -363835.80 0 -4902.3087 -410547.19 8976.0929 5688.6819 375.31742 0 -378875.89 -84832.745 724047.01 -1018090.2 0 653350.57
2300 -363855.31 0 -363855.31 0 -4893.4819 -410482.59 8977.8970 5689.5367 375.26993 0 -378898.01 -84816.035 724008.96 -1018090.9 0 653350.57
2400 -363865.13 0 -363865.13 0 -4895.2179 -410508.96 8978.6489 5689.9761 375.32747 0 -378909.08 -84812.772 723994.14 -1018090.4 0 653350.57
2500 -363868.36 0 -363868.36 0 -4894.0310 -410500.88 8978.8459 5689.3391 374.99913 0 -378911.55 -84810.894 723990.05 -1018090.7 0 653350.57
2600 -363870.68 0 -363870.68 0 -4895.0901 -410513.29 8978.7850 5689.8525 374.94772 0 -378914.26 -84813.421 723989.64 -1018090.5 0 653350.57
2700 -363871.44 0 -363871.44 0 -4894.0299 -410503.95 8979.0109 5689.7577 374.97827 0 -378915.19 -84810.089 723985.37 -1018090.5 0 653350.57
2800 -363873.17 0 -363873.17 0 -4895.6782 -410521.39 8978.8980 5690.0828 375.08553 0 -378917.24 -84814.755 723987.86 -1018090.3 0 653350.57
2900 -363874.90 0 -363874.90 0 -4898.1355 -410546.53 8978.9154 5689.0597 375.01673 0 -378917.89 -84819.463 723991.95 -1018090.4 0 653350.57
3000 -363875.80 0 -363875.80 0 -4896.9072 -410535.72 8978.9196 5689.1948 374.94777 0 -378918.86 -84818.180 723989.95 -1018090.6 0 653350.57
3100 -363877.15 0 -363877.15 0 -4895.4137 -410522.84 8978.8804 5689.5418 374.96197 0 -378920.53 -84815.238 723986.62 -1018091.9 0 653350.57
3200 -363878.66 0 -363878.66 0 -4896.9414 -410538.92 8979.0282 5689.5279 374.96731 0 -378922.19 -84815.011 723984.73 -1018091.9 0 653350.57
3272 -363878.83 0 -363878.83 0 -4897.8648 -410547.88 8978.9325 5689.3951 374.96137 0 -378922.12 -84816.968 723986.86 -1018092.0 0 653350.57
In [65]:
makeThermoPlotsFromDataFrame(minimization_thermo_pd);
In [66]:
minimization_thermo_pd[["E_long"]][2:].plot()
Out[66]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1b99839400>
In [67]:
# The total energy decreases, but intramolecular energy increases during minimization:
minimization_thermo_pd[["PotEng","E_pair"]][2:].plot()
Out[67]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1b9c56eba8>
In [68]:
# double-check: total potential energy of system minus non-bonded energy (LJ & Coulomb) 
# should correspond to intramolecular energy:
intramolecularEnergyValidation = minimization_thermo_pd["PotEng"] - minimization_thermo_pd["E_pair"]
In [69]:
intramolecularEnergyValidationDiff = (intramolecularEnergyValidation - minimization_thermo_pd["E_intramolecular"])
In [70]:
intramolecularEnergyValidationDiff.max()
Out[70]:
0.0092899999963265145
In [71]:
intramolecularEnergyValidationDiff.abs().max() / intramolecularEnergyValidation.min()
Out[71]:
7.1352862534958468e-07
In [72]:
(intramolecularEnergyValidation - minimization_thermo_pd["E_intramolecular"])[1:].plot()
# obviously "equal" (up to a tiny fraction)
Out[72]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1b9c0ba7f0>
In [73]:
# descent to steep t the first few steps, excluded
makeThermoPlotsFromDataFrame(minimization_thermo_pd.iloc[2:].copy());

NVT equilibration

In [106]:
%ls *nvtEquilibration.log
41_CTAB_on_AU_111_21x12x2_bilayer_02_nvtEquilibration.log
In [107]:
!./extract_thermo.sh *nvtEquilibration.log 02_nvtEquilibration_thermo.out
In [108]:
nvtEquilibration_thermo_file = absolute_prefix + os.sep + '02_nvtEquilibration_thermo.out'
In [109]:
!cat $nvtEquilibration_thermo_file | wc -l # count the lines in file, simple way for verification
# at 2500 steps and output every 10 steps ~ 250 steps, 251 frames, 1 header line
1002
In [110]:
nvtEquilibration_thermo_pd = pd.read_csv(nvtEquilibration_thermo_file,delim_whitespace=True)
In [111]:
nvtEquilibration_thermo_pd.set_index("Step",inplace=True)
In [112]:
nvtEquilibration_thermo_pd
Out[112]:
TotEng KinEng PotEng Temp Press Enthalpy E_bond E_angle E_dihed E_impro E_pair E_vdwl E_coul E_long E_tail Volume
Step
0 -352541.47 25886.318 -378427.79 204.92527 -2918.952200 -380354.55 7.97739 118.9759 374.96137 0 -378929.70 -84818.188 494479.45 -788590.96 0 653350.57
100 -280796.81 33114.993 -313911.80 262.15003 1992.666000 -261809.80 170.92497 1280.0829 549.44285 0 -315912.25 -89642.685 562190.12 -788459.69 0 653350.57
200 -270769.71 36365.109 -307134.82 287.87910 1142.225800 -259886.07 188.61812 1329.6683 569.61980 0 -309222.72 -93338.488 572692.23 -788576.46 0 653350.57
300 -267881.46 37079.041 -304960.50 293.53084 444.272580 -263648.23 199.06615 1318.5952 572.23374 0 -307050.39 -95399.625 577051.98 -788702.75 0 653350.57
400 -265903.64 37313.866 -303217.51 295.38980 404.875920 -262045.80 190.57540 1377.5424 610.72145 0 -305396.34 -96296.601 579666.64 -788766.38 0 653350.57
500 -264770.15 37407.605 -302177.75 296.13187 207.002240 -262797.74 177.27598 1432.3787 583.92745 0 -304371.34 -97194.261 581627.50 -788804.57 0 653350.57
600 -264348.98 37725.302 -302074.29 298.64687 204.915200 -262396.46 191.90498 1379.8082 571.44611 0 -304217.45 -97453.729 582098.81 -788862.52 0 653350.57
700 -263891.62 37709.204 -301600.82 298.51943 -29.044676 -264168.37 199.38264 1344.8841 596.37261 0 -303741.46 -98186.932 583344.77 -788899.30 0 653350.57
800 -263693.07 37669.577 -301362.65 298.20573 -40.672413 -264080.61 204.35286 1328.7876 587.45894 0 -303483.24 -98316.847 583721.88 -788888.28 0 653350.57
900 -263766.83 37711.976 -301478.81 298.54137 80.174006 -263002.90 159.53289 1396.8615 608.23865 0 -303643.44 -98351.764 583663.09 -788954.77 0 653350.57
1000 -263477.48 37884.467 -301361.95 299.90688 -151.568500 -264921.70 187.24003 1313.2100 603.57541 0 -303465.98 -98561.566 584071.85 -788976.26 0 653350.57
1100 -263719.90 37423.920 -301143.82 296.26102 -91.694352 -264593.60 183.39699 1320.0117 592.92965 0 -303240.16 -98506.333 584266.44 -789000.26 0 653350.57
1200 -263597.28 37624.391 -301221.67 297.84802 -64.449534 -264211.38 198.95474 1381.7134 593.34089 0 -303395.68 -98607.521 584221.02 -789009.18 0 653350.57
1300 -263319.03 37647.433 -300966.47 298.03043 -173.930420 -264976.32 185.38015 1334.4604 582.98327 0 -303069.29 -98645.265 584591.91 -789015.94 0 653350.57
1400 -263528.75 37585.735 -301114.48 297.54201 -131.887580 -264785.43 208.74175 1305.8240 596.98123 0 -303226.03 -98750.560 584539.47 -789014.94 0 653350.57
1500 -263661.02 37369.153 -301030.17 295.82746 -41.047229 -264052.14 210.94023 1349.3760 601.18840 0 -303191.68 -98752.496 584586.06 -789025.24 0 653350.57
1600 -263576.42 37574.612 -301151.03 297.45395 -98.555103 -264515.50 191.22945 1389.2494 607.44625 0 -303338.96 -98972.462 584666.77 -789033.27 0 653350.57
1700 -263642.84 37444.094 -301086.94 296.42072 -102.193800 -264616.59 194.84958 1346.3342 599.06474 0 -303227.19 -99007.171 584812.59 -789032.61 0 653350.57
1800 -263366.29 37845.169 -301211.46 299.59578 -154.498500 -264838.42 185.01744 1347.1790 599.68686 0 -303343.34 -99137.929 584838.56 -789043.97 0 653350.57
1900 -263540.11 37541.999 -301082.11 297.19578 -145.914070 -264930.45 202.27297 1352.7901 597.99102 0 -303235.17 -99147.266 584956.40 -789044.30 0 653350.57
2000 -263522.66 37655.787 -301178.45 298.09656 -0.252556 -263525.07 185.47487 1320.2979 595.03034 0 -303279.25 -98604.232 584375.50 -789050.52 0 653350.57
2100 -263305.13 37717.699 -301022.83 298.58668 -107.762860 -264331.94 198.75193 1392.3112 584.36427 0 -303198.26 -98950.131 584795.89 -789044.01 0 653350.57
2200 -263527.83 37554.209 -301082.04 297.29244 -100.572340 -264486.13 165.45728 1355.8058 579.72542 0 -303183.03 -98977.988 584835.83 -789040.87 0 653350.57
2300 -263742.15 37352.823 -301094.97 295.69819 -157.335730 -265241.31 179.73981 1326.9585 586.42217 0 -303188.09 -98945.903 584801.72 -789043.91 0 653350.57
2400 -263779.43 37682.182 -301461.61 298.30552 -281.760050 -266464.17 200.45882 1337.2107 596.68515 0 -303595.97 -99102.031 584553.72 -789047.66 0 653350.57
2500 -263683.99 37436.791 -301120.78 296.36291 -291.172410 -266458.41 202.81254 1360.1193 590.19938 0 -303273.91 -99296.854 585053.19 -789030.25 0 653350.57
2600 -263072.85 37787.925 -300860.77 299.14261 -298.906140 -265920.96 190.50803 1343.3646 603.96854 0 -302998.61 -99243.694 585273.34 -789028.26 0 653350.57
2700 -263468.03 37589.175 -301057.20 297.56924 -112.938960 -264544.16 176.24746 1349.9644 614.79376 0 -303198.21 -98756.928 584595.42 -789036.70 0 653350.57
2800 -263858.42 37458.519 -301316.93 296.53492 -293.314970 -266653.25 189.09428 1310.8068 602.19799 0 -303419.03 -99081.376 584704.40 -789042.05 0 653350.57
2900 -263494.71 37605.257 -301099.97 297.69655 -164.674960 -265063.81 182.46205 1351.9476 607.11154 0 -303241.49 -99179.802 584980.46 -789042.15 0 653350.57
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
97100 -263502.65 37707.116 -301209.76 298.50290 70.842987 -262827.63 191.05359 1349.9471 614.60905 0 -303365.37 -98798.435 584501.64 -789068.58 0 653350.57
97200 -263953.40 37393.230 -301346.63 296.01807 78.171734 -263208.54 183.41339 1320.1592 607.13608 0 -303457.34 -98655.899 584252.36 -789053.80 0 653350.57
97300 -263807.77 37821.839 -301629.61 299.41109 14.790298 -263666.84 181.91110 1351.0836 623.63497 0 -303786.24 -98802.227 584080.70 -789064.71 0 653350.57
97400 -264031.68 37531.730 -301563.41 297.11449 -79.597811 -264790.12 198.93731 1315.5823 614.18504 0 -303692.11 -98852.172 584215.22 -789055.16 0 653350.57
97500 -263789.70 37611.528 -301401.23 297.74620 -245.771790 -266131.52 192.11526 1371.7996 630.00891 0 -303595.15 -99457.425 584915.43 -789053.16 0 653350.57
97600 -264033.34 37510.884 -301544.23 296.94946 -229.360420 -266218.79 190.95224 1330.7553 613.87146 0 -303679.80 -99082.290 584457.31 -789054.82 0 653350.57
97700 -263709.64 37559.201 -301268.84 297.33195 -50.557978 -264191.38 200.63079 1392.9475 616.66857 0 -303479.09 -99011.161 584596.33 -789064.26 0 653350.57
97800 -264156.44 37455.575 -301612.01 296.51161 -176.962580 -265842.62 223.40239 1358.6242 606.55841 0 -303800.60 -99152.801 584422.45 -789070.24 0 653350.57
97900 -263936.53 37725.553 -301662.08 298.64886 26.702571 -263682.09 181.85293 1329.4120 611.67269 0 -303785.02 -98684.933 583978.42 -789078.51 0 653350.57
98000 -264125.17 37547.270 -301672.44 297.23751 135.593870 -262833.18 205.14994 1332.9765 614.19013 0 -303824.76 -98559.529 583798.38 -789063.62 0 653350.57
98100 -263862.26 37568.152 -301430.41 297.40281 -9.994982 -263957.50 202.88931 1351.1364 609.50722 0 -303593.95 -98905.137 584350.87 -789039.68 0 653350.57
98200 -263798.58 37600.898 -301399.48 297.66204 -126.590790 -265004.79 200.08858 1355.0205 633.33615 0 -303587.92 -98889.065 584371.12 -789069.97 0 653350.57
98300 -263610.40 37880.489 -301490.89 299.87539 28.162321 -263342.06 184.34677 1300.5267 607.59145 0 -303583.35 -98963.751 584439.42 -789059.02 0 653350.57
98400 -263607.10 37766.070 -301373.17 298.96960 -87.403957 -264439.92 194.79972 1332.7015 607.72993 0 -303508.40 -99097.428 584660.98 -789071.95 0 653350.57
98500 -263687.29 37426.293 -301113.58 296.27981 -54.873079 -264210.15 186.66144 1338.9865 611.08918 0 -303250.32 -98916.555 584714.73 -789048.50 0 653350.57
98600 -263776.63 37746.034 -301522.67 298.81099 -233.254810 -265999.19 188.17017 1335.8876 619.33498 0 -303666.06 -99374.941 584781.36 -789072.48 0 653350.57
98700 -263850.98 37754.403 -301605.38 298.87724 -121.797500 -265011.52 186.15625 1384.2293 617.14728 0 -303792.92 -98935.803 584198.96 -789056.08 0 653350.57
98800 -263710.14 37819.562 -301529.70 299.39307 -96.030522 -264625.16 191.50594 1321.8512 609.16804 0 -303652.23 -98724.989 584119.13 -789046.36 0 653350.57
98900 -263916.78 37665.239 -301582.02 298.17139 -43.304682 -264329.41 195.40419 1373.6347 609.28407 0 -303760.35 -98930.418 584233.00 -789062.93 0 653350.57
99000 -264067.59 37681.725 -301749.31 298.30190 -155.855890 -265552.65 207.94463 1286.9864 620.05725 0 -303864.30 -99033.049 584213.05 -789044.31 0 653350.57
99100 -263810.52 37496.166 -301306.69 296.83295 -322.092600 -266879.57 180.87454 1316.7644 620.32944 0 -303424.66 -99124.541 584763.30 -789063.42 0 653350.57
99200 -263927.02 37616.482 -301543.50 297.78541 45.671636 -263491.84 194.08364 1323.3538 605.63172 0 -303666.57 -98974.353 584364.47 -789056.69 0 653350.57
99300 -264000.20 37576.256 -301576.46 297.46697 17.735698 -263831.21 192.79249 1350.3039 634.62542 0 -303754.18 -98925.949 584214.13 -789042.36 0 653350.57
99400 -263878.83 37493.839 -301372.67 296.81453 -5.188889 -263928.28 162.14276 1326.2960 623.81003 0 -303484.92 -99242.028 584811.48 -789054.37 0 653350.57
99500 -263945.28 37787.238 -301732.52 299.13718 -92.964423 -264831.09 196.39208 1319.9650 586.33162 0 -303835.21 -98882.568 584103.16 -789055.80 0 653350.57
99600 -263921.26 37473.481 -301394.74 296.65337 -18.842889 -264100.80 213.90484 1386.2612 595.14864 0 -303590.05 -99057.564 584524.02 -789056.51 0 653350.57
99700 -263585.58 37600.741 -301186.32 297.66080 42.187187 -263183.60 205.73760 1337.2096 588.28073 0 -303317.55 -98997.409 584748.42 -789068.55 0 653350.57
99800 -263522.84 37561.576 -301084.41 297.35075 10.168054 -263425.95 192.03045 1338.3940 618.60681 0 -303233.44 -99114.570 584932.48 -789051.35 0 653350.57
99900 -263905.49 37443.021 -301348.51 296.41223 -73.114893 -264602.16 200.29217 1365.1563 607.68270 0 -303521.64 -99124.380 584652.23 -789049.50 0 653350.57
100000 -263879.25 37392.932 -301272.18 296.01571 -195.695670 -265743.93 198.39569 1363.6233 591.25923 0 -303425.46 -99494.266 585141.01 -789072.20 0 653350.57

1001 rows × 16 columns

In [113]:
makeThermoPlotsFromDataFrame(nvtEquilibration_thermo_pd);

NPT equilibration

In [114]:
%ls *nptEquilibration.log
41_CTAB_on_AU_111_21x12x2_bilayer_03_nptEquilibration.log
In [115]:
!./extract_thermo.sh *nptEquilibration.log 03_nptEquilibration_thermo.out
In [116]:
nptEquilibration_thermo_file = absolute_prefix + os.sep + '03_nptEquilibration_thermo.out'
In [117]:
nptEquilibration_thermo_pd = pd.read_csv(nptEquilibration_thermo_file,delim_whitespace=True)
In [118]:
nptEquilibration_thermo_pd.set_index("Step",inplace=True)
In [119]:
nptEquilibration_thermo_pd
Out[119]:
TotEng KinEng PotEng Temp Press Enthalpy E_bond E_angle E_dihed E_impro E_pair E_vdwl E_coul E_long E_tail Volume
Step
0 -263879.25 37392.932 -301272.18 296.01571 -185.402460 -265645.85 198.39569 1363.6233 591.25923 0 -303425.46 -99494.266 585141.01 -789072.20 0 653350.57
100 -263704.61 37516.284 -301220.89 296.99221 -145.116090 -265087.01 180.27196 1314.6185 605.04367 0 -303320.83 -99173.962 584918.71 -789065.58 0 653192.17
200 -263490.19 37941.826 -301432.01 300.36095 -74.301762 -264197.91 185.27422 1343.3705 618.89384 0 -303579.55 -99088.039 584572.09 -789063.61 0 653113.78
300 -263780.29 37635.831 -301416.12 297.93858 -233.936520 -266008.14 177.68944 1358.4604 623.20616 0 -303575.47 -99280.880 584753.26 -789047.86 0 652999.01
400 -263630.35 37538.717 -301169.06 297.16980 -56.808982 -264171.27 198.37089 1330.4240 620.66300 0 -303318.52 -99079.248 584814.03 -789053.30 0 652895.51
500 -263605.37 37681.599 -301286.97 298.30090 -183.653190 -265354.03 191.70493 1298.5575 611.79661 0 -303389.03 -99058.319 584723.57 -789054.28 0 652874.55
600 -263723.18 37775.581 -301498.76 299.04489 -152.129070 -265171.54 184.70303 1363.6115 599.68747 0 -303646.76 -98887.643 584280.45 -789039.57 0 652812.79
700 -264032.14 37590.247 -301622.39 297.57772 -84.774937 -264839.21 196.66229 1334.5266 584.02648 0 -303737.61 -98786.182 584115.75 -789067.17 0 652777.15
800 -263799.61 37421.683 -301221.29 296.24331 -159.309620 -265316.05 181.62709 1386.3291 618.01179 0 -303407.26 -98903.637 584573.99 -789077.61 0 652689.78
900 -263769.34 37549.128 -301318.47 297.25221 -175.172370 -265436.40 182.41557 1357.4513 605.76605 0 -303464.10 -99019.951 584615.44 -789059.59 0 652544.19
1000 -263972.46 37476.076 -301448.54 296.67391 -25.766700 -264217.64 185.24997 1336.9823 583.81612 0 -303554.59 -98786.998 584284.48 -789052.07 0 652455.16
1100 -264294.98 37692.832 -301987.81 298.38983 -210.165250 -266294.52 179.19113 1336.5248 617.72558 0 -304121.25 -99047.332 583982.71 -789056.63 0 652370.73
1200 -263584.34 37818.581 -301402.92 299.38530 -139.099420 -264907.56 195.48901 1344.0313 617.33452 0 -303559.78 -99152.962 584632.72 -789039.53 0 652274.87
1300 -263779.39 37694.544 -301473.93 298.40338 -346.706710 -267076.90 203.00633 1364.6221 611.31498 0 -303652.87 -99206.154 584618.73 -789065.45 0 652152.26
1400 -263792.61 37788.827 -301581.43 299.14975 -90.487318 -264653.09 180.30340 1336.6452 611.83677 0 -303710.22 -98957.535 584323.16 -789075.84 0 652045.93
1500 -264043.99 37652.848 -301696.84 298.07330 -29.615160 -264325.58 179.09953 1374.9411 604.09836 0 -303854.98 -98948.824 584161.32 -789067.48 0 651968.92
1600 -264126.65 37684.788 -301811.44 298.32615 -78.653810 -264874.43 180.74649 1291.7220 597.53753 0 -303881.44 -98862.070 584039.61 -789058.99 0 651898.53
1700 -264113.79 37660.924 -301774.72 298.13723 -55.471350 -264641.16 205.17661 1317.4074 598.80735 0 -303896.11 -99106.394 584275.18 -789064.90 0 651884.48
1800 -264255.70 37520.115 -301775.82 297.02253 -167.080590 -265844.03 206.15866 1329.9045 595.11394 0 -303907.00 -99094.067 584249.17 -789062.10 0 651835.33
1900 -264037.16 37603.868 -301641.03 297.68556 -164.195920 -265597.82 175.53100 1322.8482 613.56397 0 -303752.97 -99009.425 584323.04 -789066.59 0 651734.63
2000 -264133.30 37567.731 -301701.03 297.39948 -56.594492 -264671.18 175.39320 1315.6885 601.47585 0 -303793.59 -99027.044 584308.10 -789074.64 0 651680.59
2100 -263970.95 37479.515 -301450.46 296.70114 -23.522459 -264194.47 175.97468 1348.2436 622.07357 0 -303596.75 -98835.969 584294.75 -789055.53 0 651590.67
2200 -264072.57 37488.555 -301561.12 296.77270 -245.614050 -266406.04 189.66556 1293.8347 610.51594 0 -303655.14 -99190.051 584598.38 -789063.47 0 651439.97
2300 -263934.77 37475.173 -301409.94 296.66676 -89.486017 -264784.75 170.73724 1328.0518 595.09063 0 -303503.82 -99019.944 584563.05 -789046.92 0 651296.29
2400 -263832.18 37572.633 -301404.81 297.43828 -337.546850 -267037.60 205.81106 1332.0391 613.32826 0 -303555.99 -99495.256 585003.49 -789064.23 0 651141.28
2500 -264005.30 37552.592 -301557.90 297.27963 -219.827990 -266092.46 195.20674 1336.1385 612.67794 0 -303701.92 -99139.174 584499.88 -789062.63 0 651022.26
2600 -263969.35 37645.530 -301614.88 298.01536 -209.271060 -265955.80 201.09970 1357.8026 603.00158 0 -303776.78 -99252.728 584556.60 -789080.65 0 650866.99
2700 -263919.59 37451.346 -301370.93 296.47814 -154.255800 -265383.61 202.13044 1370.9327 628.21927 0 -303572.22 -98996.530 584487.95 -789063.63 0 650774.43
2800 -263814.58 37670.953 -301485.53 298.21662 -187.653070 -265595.23 201.87031 1353.9607 592.92272 0 -303634.28 -99086.893 584512.01 -789059.40 0 650649.71
2900 -264120.85 37477.509 -301598.36 296.68525 -79.885022 -264878.79 203.01157 1376.9646 608.50743 0 -303786.85 -98953.199 584231.99 -789065.63 0 650564.65
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
97100 -265543.10 37896.867 -303439.97 300.00504 104.226110 -264587.66 207.81106 1308.5833 594.70369 0 -305551.07 -98153.999 581684.11 -789081.18 0 628566.97
97200 -265744.50 37442.603 -303187.10 296.40892 -96.562119 -266629.67 204.65556 1378.9643 594.61766 0 -305365.34 -98734.473 582460.06 -789090.93 0 628552.91
97300 -265875.61 37674.274 -303549.88 298.24292 -61.801524 -266442.09 177.35951 1319.3401 607.40453 0 -305653.98 -98608.283 582019.69 -789065.39 0 628511.59
97400 -265983.30 37513.871 -303497.17 296.97311 -55.071280 -266488.09 200.53926 1300.3425 632.47391 0 -305630.52 -98572.282 582027.55 -789085.79 0 628515.04
97500 -265885.82 37441.728 -303327.55 296.40200 -122.524540 -267008.77 184.21146 1328.3281 589.47402 0 -305429.56 -98513.762 582161.64 -789077.44 0 628436.55
97600 -266011.32 37366.637 -303377.96 295.80755 -134.829520 -267246.94 184.48564 1353.2760 612.10826 0 -305527.83 -98445.706 581987.18 -789069.31 0 628381.07
97700 -265593.10 37705.392 -303298.49 298.48925 -135.018070 -266830.26 196.10443 1311.2037 616.23616 0 -305422.04 -98767.251 582432.81 -789087.59 0 628285.50
97800 -265321.87 37522.189 -302844.06 297.03896 46.443924 -264896.33 197.28639 1346.1136 620.78205 0 -305008.24 -98644.690 582706.27 -789069.82 0 628250.56
97900 -265014.35 37958.379 -302972.73 300.49199 78.095487 -264298.73 196.07253 1399.0573 615.27935 0 -305183.14 -98484.582 582365.03 -789063.59 0 628315.67
98000 -265332.28 37737.157 -303069.44 298.74072 -96.158199 -266213.40 195.05657 1354.4192 609.47883 0 -305228.39 -98851.577 582699.48 -789076.29 0 628311.67
98100 -265619.65 37628.998 -303248.65 297.88449 -218.141120 -267618.44 198.09162 1319.3374 631.23231 0 -305397.31 -98914.402 582585.09 -789068.00 0 628280.31
98200 -265513.66 37612.096 -303125.75 297.75069 -117.440450 -266589.71 179.07190 1377.3625 608.32754 0 -305290.51 -98553.775 582346.98 -789083.72 0 628258.88
98300 -265413.13 37648.455 -303061.58 298.03852 -247.030230 -267676.10 198.42758 1427.0209 637.20935 0 -305324.24 -99031.508 582781.29 -789074.02 0 628135.78
98400 -265256.17 37790.267 -303046.44 299.16116 -126.278630 -266412.81 195.34691 1391.6244 613.44175 0 -305246.85 -98615.919 582426.14 -789057.07 0 628045.53
98500 -265416.03 37674.743 -303090.77 298.24663 71.032800 -264765.44 200.86479 1396.4696 618.34283 0 -305306.45 -98561.043 582339.81 -789085.22 0 628016.59
98600 -265977.24 37580.340 -303557.58 297.49930 208.285010 -264069.40 181.71727 1310.7394 610.96198 0 -305661.00 -98433.610 581872.01 -789099.40 0 628071.29
98700 -265883.48 37453.662 -303337.14 296.49647 1.231729 -265872.20 190.12067 1370.5156 642.87566 0 -305540.66 -98592.496 582141.27 -789089.43 0 628110.88
98800 -265867.78 37577.137 -303444.92 297.47394 56.861693 -265346.90 199.01080 1331.3026 628.81873 0 -305604.05 -98560.842 582050.57 -789093.78 0 628127.34
98900 -265496.61 37892.941 -303389.55 299.97396 102.299280 -264559.39 191.79214 1368.7296 631.41386 0 -305581.49 -98686.849 582175.47 -789070.10 0 628189.98
99000 -265517.54 37883.012 -303400.55 299.89536 -19.390457 -265695.21 176.14942 1307.9065 643.53943 0 -305528.15 -98683.356 582242.46 -789087.25 0 628262.36
99100 -265503.12 37618.058 -303121.18 297.79789 92.584285 -264654.80 190.61378 1342.4459 653.80806 0 -305308.05 -98298.185 582078.51 -789088.37 0 628269.25
99200 -265784.81 37525.861 -303310.67 297.06802 -9.503912 -265871.89 185.22759 1347.7817 624.87145 0 -305468.55 -98734.784 582357.61 -789091.38 0 628266.90
99300 -265369.17 37596.075 -302965.25 297.62386 81.949749 -264618.26 190.41732 1364.7511 632.61329 0 -305153.03 -98461.327 582383.53 -789075.24 0 628296.23
99400 -265044.46 37957.653 -303002.11 300.48624 137.149820 -263787.70 175.11676 1355.3943 633.16074 0 -305165.79 -98369.762 582280.40 -789076.43 0 628322.71
99500 -265363.35 37924.038 -303287.39 300.22013 -1.185586 -265374.21 184.69305 1322.8084 628.30540 0 -305423.19 -98421.738 582090.92 -789092.38 0 628352.43
99600 -265266.65 37826.403 -303093.05 299.44722 -46.420767 -265692.05 178.78190 1348.5517 625.00287 0 -305245.39 -98394.530 582224.51 -789075.37 0 628369.05
99700 -265611.22 37551.949 -303163.17 297.27455 -25.395548 -265843.95 181.37656 1329.6313 618.39413 0 -305292.57 -98352.778 582147.87 -789087.67 0 628363.71
99800 -265550.75 37574.484 -303125.23 297.45294 126.926610 -264387.51 194.43993 1358.3964 625.10289 0 -305303.17 -98145.022 581909.20 -789067.35 0 628405.75
99900 -265442.18 37564.613 -303006.80 297.37480 145.181390 -264111.54 193.96456 1338.1776 615.69053 0 -305154.63 -98254.383 582170.48 -789070.72 0 628454.88
100000 -265255.19 37889.409 -303144.60 299.94600 -42.324388 -265643.14 196.44234 1340.6013 630.73466 0 -305312.38 -98429.549 582194.60 -789077.43 0 628497.72

1001 rows × 16 columns

In [120]:
nptEquilibration_thermo_consecutive_pd = nptEquilibration_thermo_pd.copy()
In [121]:
nvtEquilibration_thermo_pd.index.max()
Out[121]:
100000
In [122]:
nptEquilibration_thermo_consecutive_pd.index = nptEquilibration_thermo_pd.index + nvtEquilibration_thermo_pd.index.max()
In [123]:
equilibration_pd = pd.concat([nvtEquilibration_thermo_pd,nptEquilibration_thermo_consecutive_pd])
In [124]:
equilibration_pd
Out[124]:
E_angle E_bond E_coul E_dihed E_impro E_intramolecular E_long E_pair E_tail E_vdwl Enthalpy KinEng PotEng Press Temp TotEng Volume
Step
0 118.9759 7.97739 494479.45 374.96137 0 501.91466 -788590.96 -378929.70 0 -84818.188 -380354.55 25886.318 -378427.79 -2918.952200 204.92527 -352541.47 653350.57
100 1280.0829 170.92497 562190.12 549.44285 0 2000.45072 -788459.69 -315912.25 0 -89642.685 -261809.80 33114.993 -313911.80 1992.666000 262.15003 -280796.81 653350.57
200 1329.6683 188.61812 572692.23 569.61980 0 2087.90622 -788576.46 -309222.72 0 -93338.488 -259886.07 36365.109 -307134.82 1142.225800 287.87910 -270769.71 653350.57
300 1318.5952 199.06615 577051.98 572.23374 0 2089.89509 -788702.75 -307050.39 0 -95399.625 -263648.23 37079.041 -304960.50 444.272580 293.53084 -267881.46 653350.57
400 1377.5424 190.57540 579666.64 610.72145 0 2178.83925 -788766.38 -305396.34 0 -96296.601 -262045.80 37313.866 -303217.51 404.875920 295.38980 -265903.64 653350.57
500 1432.3787 177.27598 581627.50 583.92745 0 2193.58213 -788804.57 -304371.34 0 -97194.261 -262797.74 37407.605 -302177.75 207.002240 296.13187 -264770.15 653350.57
600 1379.8082 191.90498 582098.81 571.44611 0 2143.15929 -788862.52 -304217.45 0 -97453.729 -262396.46 37725.302 -302074.29 204.915200 298.64687 -264348.98 653350.57
700 1344.8841 199.38264 583344.77 596.37261 0 2140.63935 -788899.30 -303741.46 0 -98186.932 -264168.37 37709.204 -301600.82 -29.044676 298.51943 -263891.62 653350.57
800 1328.7876 204.35286 583721.88 587.45894 0 2120.59940 -788888.28 -303483.24 0 -98316.847 -264080.61 37669.577 -301362.65 -40.672413 298.20573 -263693.07 653350.57
900 1396.8615 159.53289 583663.09 608.23865 0 2164.63304 -788954.77 -303643.44 0 -98351.764 -263002.90 37711.976 -301478.81 80.174006 298.54137 -263766.83 653350.57
1000 1313.2100 187.24003 584071.85 603.57541 0 2104.02544 -788976.26 -303465.98 0 -98561.566 -264921.70 37884.467 -301361.95 -151.568500 299.90688 -263477.48 653350.57
1100 1320.0117 183.39699 584266.44 592.92965 0 2096.33834 -789000.26 -303240.16 0 -98506.333 -264593.60 37423.920 -301143.82 -91.694352 296.26102 -263719.90 653350.57
1200 1381.7134 198.95474 584221.02 593.34089 0 2174.00903 -789009.18 -303395.68 0 -98607.521 -264211.38 37624.391 -301221.67 -64.449534 297.84802 -263597.28 653350.57
1300 1334.4604 185.38015 584591.91 582.98327 0 2102.82382 -789015.94 -303069.29 0 -98645.265 -264976.32 37647.433 -300966.47 -173.930420 298.03043 -263319.03 653350.57
1400 1305.8240 208.74175 584539.47 596.98123 0 2111.54698 -789014.94 -303226.03 0 -98750.560 -264785.43 37585.735 -301114.48 -131.887580 297.54201 -263528.75 653350.57
1500 1349.3760 210.94023 584586.06 601.18840 0 2161.50463 -789025.24 -303191.68 0 -98752.496 -264052.14 37369.153 -301030.17 -41.047229 295.82746 -263661.02 653350.57
1600 1389.2494 191.22945 584666.77 607.44625 0 2187.92510 -789033.27 -303338.96 0 -98972.462 -264515.50 37574.612 -301151.03 -98.555103 297.45395 -263576.42 653350.57
1700 1346.3342 194.84958 584812.59 599.06474 0 2140.24852 -789032.61 -303227.19 0 -99007.171 -264616.59 37444.094 -301086.94 -102.193800 296.42072 -263642.84 653350.57
1800 1347.1790 185.01744 584838.56 599.68686 0 2131.88330 -789043.97 -303343.34 0 -99137.929 -264838.42 37845.169 -301211.46 -154.498500 299.59578 -263366.29 653350.57
1900 1352.7901 202.27297 584956.40 597.99102 0 2153.05409 -789044.30 -303235.17 0 -99147.266 -264930.45 37541.999 -301082.11 -145.914070 297.19578 -263540.11 653350.57
2000 1320.2979 185.47487 584375.50 595.03034 0 2100.80311 -789050.52 -303279.25 0 -98604.232 -263525.07 37655.787 -301178.45 -0.252556 298.09656 -263522.66 653350.57
2100 1392.3112 198.75193 584795.89 584.36427 0 2175.42740 -789044.01 -303198.26 0 -98950.131 -264331.94 37717.699 -301022.83 -107.762860 298.58668 -263305.13 653350.57
2200 1355.8058 165.45728 584835.83 579.72542 0 2100.98850 -789040.87 -303183.03 0 -98977.988 -264486.13 37554.209 -301082.04 -100.572340 297.29244 -263527.83 653350.57
2300 1326.9585 179.73981 584801.72 586.42217 0 2093.12048 -789043.91 -303188.09 0 -98945.903 -265241.31 37352.823 -301094.97 -157.335730 295.69819 -263742.15 653350.57
2400 1337.2107 200.45882 584553.72 596.68515 0 2134.35467 -789047.66 -303595.97 0 -99102.031 -266464.17 37682.182 -301461.61 -281.760050 298.30552 -263779.43 653350.57
2500 1360.1193 202.81254 585053.19 590.19938 0 2153.13122 -789030.25 -303273.91 0 -99296.854 -266458.41 37436.791 -301120.78 -291.172410 296.36291 -263683.99 653350.57
2600 1343.3646 190.50803 585273.34 603.96854 0 2137.84117 -789028.26 -302998.61 0 -99243.694 -265920.96 37787.925 -300860.77 -298.906140 299.14261 -263072.85 653350.57
2700 1349.9644 176.24746 584595.42 614.79376 0 2141.00562 -789036.70 -303198.21 0 -98756.928 -264544.16 37589.175 -301057.20 -112.938960 297.56924 -263468.03 653350.57
2800 1310.8068 189.09428 584704.40 602.19799 0 2102.09907 -789042.05 -303419.03 0 -99081.376 -266653.25 37458.519 -301316.93 -293.314970 296.53492 -263858.42 653350.57
2900 1351.9476 182.46205 584980.46 607.11154 0 2141.52119 -789042.15 -303241.49 0 -99179.802 -265063.81 37605.257 -301099.97 -164.674960 297.69655 -263494.71 653350.57
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
197100 1308.5833 207.81106 581684.11 594.70369 0 NaN -789081.18 -305551.07 0 -98153.999 -264587.66 37896.867 -303439.97 104.226110 300.00504 -265543.10 628566.97
197200 1378.9643 204.65556 582460.06 594.61766 0 NaN -789090.93 -305365.34 0 -98734.473 -266629.67 37442.603 -303187.10 -96.562119 296.40892 -265744.50 628552.91
197300 1319.3401 177.35951 582019.69 607.40453 0 NaN -789065.39 -305653.98 0 -98608.283 -266442.09 37674.274 -303549.88 -61.801524 298.24292 -265875.61 628511.59
197400 1300.3425 200.53926 582027.55 632.47391 0 NaN -789085.79 -305630.52 0 -98572.282 -266488.09 37513.871 -303497.17 -55.071280 296.97311 -265983.30 628515.04
197500 1328.3281 184.21146 582161.64 589.47402 0 NaN -789077.44 -305429.56 0 -98513.762 -267008.77 37441.728 -303327.55 -122.524540 296.40200 -265885.82 628436.55
197600 1353.2760 184.48564 581987.18 612.10826 0 NaN -789069.31 -305527.83 0 -98445.706 -267246.94 37366.637 -303377.96 -134.829520 295.80755 -266011.32 628381.07
197700 1311.2037 196.10443 582432.81 616.23616 0 NaN -789087.59 -305422.04 0 -98767.251 -266830.26 37705.392 -303298.49 -135.018070 298.48925 -265593.10 628285.50
197800 1346.1136 197.28639 582706.27 620.78205 0 NaN -789069.82 -305008.24 0 -98644.690 -264896.33 37522.189 -302844.06 46.443924 297.03896 -265321.87 628250.56
197900 1399.0573 196.07253 582365.03 615.27935 0 NaN -789063.59 -305183.14 0 -98484.582 -264298.73 37958.379 -302972.73 78.095487 300.49199 -265014.35 628315.67
198000 1354.4192 195.05657 582699.48 609.47883 0 NaN -789076.29 -305228.39 0 -98851.577 -266213.40 37737.157 -303069.44 -96.158199 298.74072 -265332.28 628311.67
198100 1319.3374 198.09162 582585.09 631.23231 0 NaN -789068.00 -305397.31 0 -98914.402 -267618.44 37628.998 -303248.65 -218.141120 297.88449 -265619.65 628280.31
198200 1377.3625 179.07190 582346.98 608.32754 0 NaN -789083.72 -305290.51 0 -98553.775 -266589.71 37612.096 -303125.75 -117.440450 297.75069 -265513.66 628258.88
198300 1427.0209 198.42758 582781.29 637.20935 0 NaN -789074.02 -305324.24 0 -99031.508 -267676.10 37648.455 -303061.58 -247.030230 298.03852 -265413.13 628135.78
198400 1391.6244 195.34691 582426.14 613.44175 0 NaN -789057.07 -305246.85 0 -98615.919 -266412.81 37790.267 -303046.44 -126.278630 299.16116 -265256.17 628045.53
198500 1396.4696 200.86479 582339.81 618.34283 0 NaN -789085.22 -305306.45 0 -98561.043 -264765.44 37674.743 -303090.77 71.032800 298.24663 -265416.03 628016.59
198600 1310.7394 181.71727 581872.01 610.96198 0 NaN -789099.40 -305661.00 0 -98433.610 -264069.40 37580.340 -303557.58 208.285010 297.49930 -265977.24 628071.29
198700 1370.5156 190.12067 582141.27 642.87566 0 NaN -789089.43 -305540.66 0 -98592.496 -265872.20 37453.662 -303337.14 1.231729 296.49647 -265883.48 628110.88
198800 1331.3026 199.01080 582050.57 628.81873 0 NaN -789093.78 -305604.05 0 -98560.842 -265346.90 37577.137 -303444.92 56.861693 297.47394 -265867.78 628127.34
198900 1368.7296 191.79214 582175.47 631.41386 0 NaN -789070.10 -305581.49 0 -98686.849 -264559.39 37892.941 -303389.55 102.299280 299.97396 -265496.61 628189.98
199000 1307.9065 176.14942 582242.46 643.53943 0 NaN -789087.25 -305528.15 0 -98683.356 -265695.21 37883.012 -303400.55 -19.390457 299.89536 -265517.54 628262.36
199100 1342.4459 190.61378 582078.51 653.80806 0 NaN -789088.37 -305308.05 0 -98298.185 -264654.80 37618.058 -303121.18 92.584285 297.79789 -265503.12 628269.25
199200 1347.7817 185.22759 582357.61 624.87145 0 NaN -789091.38 -305468.55 0 -98734.784 -265871.89 37525.861 -303310.67 -9.503912 297.06802 -265784.81 628266.90
199300 1364.7511 190.41732 582383.53 632.61329 0 NaN -789075.24 -305153.03 0 -98461.327 -264618.26 37596.075 -302965.25 81.949749 297.62386 -265369.17 628296.23
199400 1355.3943 175.11676 582280.40 633.16074 0 NaN -789076.43 -305165.79 0 -98369.762 -263787.70 37957.653 -303002.11 137.149820 300.48624 -265044.46 628322.71
199500 1322.8084 184.69305 582090.92 628.30540 0 NaN -789092.38 -305423.19 0 -98421.738 -265374.21 37924.038 -303287.39 -1.185586 300.22013 -265363.35 628352.43
199600 1348.5517 178.78190 582224.51 625.00287 0 NaN -789075.37 -305245.39 0 -98394.530 -265692.05 37826.403 -303093.05 -46.420767 299.44722 -265266.65 628369.05
199700 1329.6313 181.37656 582147.87 618.39413 0 NaN -789087.67 -305292.57 0 -98352.778 -265843.95 37551.949 -303163.17 -25.395548 297.27455 -265611.22 628363.71
199800 1358.3964 194.43993 581909.20 625.10289 0 NaN -789067.35 -305303.17 0 -98145.022 -264387.51 37574.484 -303125.23 126.926610 297.45294 -265550.75 628405.75
199900 1338.1776 193.96456 582170.48 615.69053 0 NaN -789070.72 -305154.63 0 -98254.383 -264111.54 37564.613 -303006.80 145.181390 297.37480 -265442.18 628454.88
200000 1340.6013 196.44234 582194.60 630.73466 0 NaN -789077.43 -305312.38 0 -98429.549 -265643.14 37889.409 -303144.60 -42.324388 299.94600 -265255.19 628497.72

2002 rows × 17 columns

In [125]:
makeThermoPlotsFromDataFrame(nptEquilibration_thermo_pd);
In [126]:
# show nvt and npt equilibration consecutively
makeThermoPlotsFromDataFrame(equilibration_pd);
In [127]:
makeRollingAverageThermoPlotsFromDataFrame(equilibration_pd,window=10);

1 ns NPT

In [128]:
%ls -lh *_1ns_npt_01_nptProduction.log
ls: cannot access *_1ns_npt_01_nptProduction.log: No such file or directory
In [61]:
!./extract_thermo.sh *_1ns_npt_01_nptProduction.log 05_1ns_nptProduction_thermo.out
In [62]:
nptProduction_1ns_thermo_file = absolute_prefix + os.sep + '05_1ns_nptProduction_thermo.out'
nptProduction_1ns_thermo_pd = pd.read_csv(nptProduction_1ns_thermo_file,delim_whitespace=True)
nptProduction_1ns_thermo_pd.set_index("Step",inplace=True)
In [63]:
makeThermoPlotsFromDataFrame(nptProduction_1ns_thermo_pd.iloc[::100].copy()); # only every 100th data point
In [64]:
makeRollingAverageThermoPlotsFromDataFrame(nptProduction_1ns_thermo_pd,window=5000);
In [65]:
makeRollingAverageThermoPlotsFromDataFrame(nptProduction_1ns_thermo_pd,window=25000);
In [66]:
nptProduction_1ns_thermo_pd[["PotEng","E_pair"]].rolling(window=5000,center=True).mean().plot()
Out[66]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa85676ea90>

1 ns NVE

In [67]:
%ls -lh *_1ns_nve_01_nveProduction.log
-rw-r--r-- 1 fr_jh1130 fr_fr 11M Jul 27 06:04 1_CTAB_on_AU_111_21x12x2_1ns_nve_01_nveProduction.log
In [68]:
!./extract_thermo.sh *_1ns_nve_01_nveProduction.log 05_1ns_nveProduction_thermo.out
In [69]:
nveProduction_1ns_thermo_file = absolute_prefix + os.sep + '05_1ns_nveProduction_thermo.out'
nveProduction_1ns_thermo_pd = pd.read_csv(nveProduction_1ns_thermo_file,delim_whitespace=True)
nveProduction_1ns_thermo_pd.set_index("Step",inplace=True)
In [70]:
makeThermoPlotsFromDataFrame(nveProduction_1ns_thermo_pd.iloc[::100].copy()); # only every 100th data point
In [71]:
makeRollingAverageThermoPlotsFromDataFrame(nveProduction_1ns_thermo_pd,window=5000);
In [72]:
makeRollingAverageThermoPlotsFromDataFrame(nveProduction_1ns_thermo_pd,window=25000);
In [73]:
nveProduction_1ns_thermo_pd[["PotEng","E_pair"]].rolling(window=5000,center=True).mean().plot()
Out[73]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa858fff7b8>
In [74]:
ax = nveProduction_1ns_thermo_pd[["TotEng"]].rolling(window=5000,center=True).mean().plot()
nptProduction_1ns_thermo_pd[["TotEng"]].rolling(window=5000,center=True).mean().plot(ax=ax)
Out[74]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa85623fcf8>
In [75]:
nveProduction_1ns_thermo_pd[["TotEng"]].rolling(window=5000,center=True).mean().plot()
Out[75]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa85675eda0>
In [76]:
nptProduction_1ns_thermo_pd[["TotEng"]].rolling(window=5000,center=True).mean().plot()
Out[76]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa85623f160>

Trajectory visualization with ASE and ParmEd

In [77]:
%ls *.lammps
1_CTAB_on_AU_111_21x12x2_1ns_npt_final_config.lammps
1_CTAB_on_AU_111_21x12x2_1ns_nve_final_config.lammps
1_CTAB_on_AU_111_21x12x2_minimized.lammps
1_CTAB_on_AU_111_21x12x2_nptEquilibrated.lammps
1_CTAB_on_AU_111_21x12x2_nvtEquilibrated.lammps
log.lammps
In [78]:
%ls -lh *.nc
-rw-r--r-- 1 fr_jh1130 fr_fr 5.5G Jul 27 03:16 1_CTAB_on_AU_111_21x12x2_1ns_npt_nptProduction.nc
-rw-r--r-- 1 fr_jh1130 fr_fr 8.2G Jul 27 06:04 1_CTAB_on_AU_111_21x12x2_1ns_nve_nveProduction.nc
-rw-r--r-- 1 fr_jh1130 fr_fr  29M Jul 26 19:30 1_CTAB_on_AU_111_21x12x2_nptEquilibration.nc
-rw-r--r-- 1 fr_jh1130 fr_fr  29M Jul 26 19:26 1_CTAB_on_AU_111_21x12x2_nvtEquilibration.nc
In [79]:
# File names
prefix = '1_CTAB_on_AU_111_21x12x2_'

pdb_file_initial_config = absolute_prefix + os.sep + prefix + 'psfgen_ctrl.pdb'

lmp_files = { 'initial':         prefix + 'psfgen_CLA2BR.data',
              'minimized':       prefix + 'minimized.lammps',
              'nvtEquilibrated': prefix + 'nvtEquilibrated.lammps',
              'nptEquilibrated': prefix + 'nptEquilibrated.lammps',
              #'npt10ps':         prefix + '10ps_npt_final_config.lammps',
              #'npt100ps':        prefix + '100ps_npt_final_config.lammps',
              'npt1ns':          prefix + '1ns_npt_final_config.lammps',
              'nve1ns':          prefix + '1ns_nve_final_config.lammps'
            }

# lmp_dumps = { 
#              'nvtEquilibration': prefix + 'nvtEquilibration.dump',
#              #'nptEquilibration': prefix + 'nptEquilibration.dump',
#              #'npt10ps':          prefix + '10ps_npt_nptProduction.dump',
#              #'npt100ps':         prefix + '100ps_npt_nptProduction.dump',
#              #'npt1ns':           prefix + '1ns_npt_with_restarts_nptProduction.dump',
#              #'nve1ns':           prefix + '1ns_nve_with_restarts_nveProduction.dump' 
#             }
lmp_netcdf = { 
              'nvtEquilibration': prefix + 'nvtEquilibration.nc',
              'nptEquilibration': prefix + 'nptEquilibration.nc',
              'npt1ns':           prefix + '1ns_npt_nptProduction.nc',
              'nve1ns':           prefix + '1ns_nve_nveProduction.nc' 
             }
In [80]:
# construct a dictionary-like atom type-> atom number array
t2n, t2n_array = inferTypes2NumbersFromPdbAndLmp(pdb_file_initial_config, lmp_files['initial'])
PDB contains following residue types ['AUM' 'CTA' 'HOH' 'ION']
System contains 11 atom types.
In [81]:
t2n
Out[81]:
{1: 1, 2: 1, 3: 1, 4: 6, 5: 6, 6: 6, 7: 7, 8: 79, 9: 1, 10: 8, 11: 35}
In [82]:
t2n_array
Out[82]:
array([ 0,  1,  1,  1,  6,  6,  6,  7, 79,  1,  8, 35], dtype=uint64)
In [83]:
np.array(ase.data.chemical_symbols)[t2n_array] # double-check against LAMMPS data file
Out[83]:
array(['X', 'H', 'H', 'H', 'C', 'C', 'C', 'N', 'Au', 'H', 'O', 'Br'],
      dtype='<U2')

Initial configuration from .pdb

In [84]:
# create atom selections for later post-processing:
struct_pdb_pmd     = pmd.read_PDB(pdb_file_initial_config)
water              = [ a for a in struct_pdb_pmd.atoms if a.residue.name == 'HOH' ] # water
water_indices      = [ a.number - 1 for a in water ] # to remove water atoms later
surface            = [ a for a in struct_pdb_pmd.atoms if a.residue.name == 'AUM' ] # gold surface
surface_indices    = [ a.number - 1 for a in surface ] 
surfactant         = [ a for a in struct_pdb_pmd.atoms if a.residue.name == 'CTA' ]
surfactant_indices = [ a.number - 1 for a in surfactant ] 
ions               = [ a for a in struct_pdb_pmd.atoms if a.residue.name == 'ION' ] # sodium counterions
ion_indices        = [ a.number - 1 for a in ions ] 
In [85]:
struct_pdb_pmd
Out[85]:
<Structure 48702 atoms; 18231 residues; 30471 bonds; NOT parametrized>
In [86]:
surface_indices[0] # double check index begins at 0
Out[86]:
0
In [87]:
len(surfactant_indices)
Out[87]:
62
In [88]:
ions
Out[88]:
[<Atom CLA [16556]; In ION 7515>]
In [89]:
ions[0].atomic_number # wrong: inferred C (carbon)
Out[89]:
6
In [90]:
nv.show_parmed(struct_pdb_pmd)

LAMMPS data files

In [94]:
# read frames of interest
lmp_frames = {}
for k,f in lmp_files.items():
    lmp_frames[k] = read(lmp_files[k],format='lammps-data')
    lmp_frames[k].set_atomic_numbers(
        t2n_array[lmp_frames[k].get_atomic_numbers() ] )
In [1]:
lmp_frames['minimized']
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-08a4878d1c3b> in <module>()
----> 1 lmp_frames

NameError: name 'lmp_frames' is not defined
In [96]:
# lmp_views = []
# for k, f in lmp_frames.items():
#     lmp_views.append( nv.show_ase(f) )
#     lmp_views[-1]._set_sync_camera()
#     lmp_views[-1]._remote_call("setSize", target="Widget", args=["250px", "250px"])
#     lmp_views[-1].center()
#     lmp_views[-1].render_image()
In [97]:
# vbox = ipywidgets.VBox(lmp_views)
In [98]:
# vbox
In [99]:
# frames stripped of all oxygen & hydrogen:
# apparently, order of atoms in dump trajectories is fixed, but 
# in LAMMPS data files not, hence the indexing does not work based upon previously derived sets
lmp_naked_frames = {}
for k,f in lmp_frames.items():
    #lmp_naked_frames[k] = []
    #or f in lmp_frames[k]:
    #g = f.copy()
    g = f[ (np.array(f.get_chemical_symbols()) != 'O') & (np.array(f.get_chemical_symbols()) != 'H')]
    #del g[water_indices]
    lmp_naked_frames[k] = g.copy()
In [100]:
lmp_naked_frames["nvtEquilibrated"]
Out[100]:
Atoms(symbols='C19Au3024BrN', pbc=True, cell=[64.107506, 64.107506, 128.215012], angles=..., bonds=..., id=..., masses=..., mol-id=..., momenta=..., travel=..., type=...)
In [101]:
lmp_frames["nvtEquilibrated"]
Out[101]:
Atoms(symbols='C19H30452Au3024BrNO15205', pbc=True, cell=[64.107506, 64.107506, 128.215012], angles=..., bonds=..., id=..., masses=..., mol-id=..., momenta=..., travel=..., type=...)
In [102]:
lmp_views = []
for k, f in lmp_naked_frames.items():
    lmp_views.append( nv.show_ase(f) )
    lmp_views[-1]._set_sync_camera()
    lmp_views[-1]._remote_call("setSize", target="Widget", args=["250px", "250px"])
    lmp_views[-1].center()
    lmp_views[-1].render_image()
In [103]:
vbox = ipywidgets.VBox(lmp_views)
In [104]:
vbox

LAMMPS NETCDF trajectories

In [106]:
# lmp_trajectrories = {}
# for k,t in lmp_dumps.items():
#     lmp_trajectrories[k] = read(t, index=':',format='lammps-dump')
#     for f in lmp_trajectrories[k]:      
#         f.set_atomic_numbers(
#             t2n_array[f.get_atomic_numbers() ] )
#         f.center()
In [107]:
# lmp_netcdf['nvtEquilibration']
In [108]:
# # trial
# traj = NetCDFTrajectory(lmp_netcdf['nvtEquilibration'], 'r',
#                         types_to_numbers=list(t2n_array),
#                         keep_open=True)
# # workaround to center trajecory in box
# # (ASE places the box origin at 0,0,0 by default)
# frames = []
# for i,frame in enumerate(traj):
#     frame.center() 
#     frames.append( frame )
# traj.close()
In [109]:
lmp_trajectrories = {}
for k,t in lmp_netcdf.items():
    #lmp_trajectrories[k] = read(t, index=':',format='lammps-dump')
    tmp_traj = NetCDFTrajectory(lmp_netcdf[k], 'r',
                        types_to_numbers=list(t2n_array),
                        keep_open=True)
    # workaround to center trajecory in box
    # (ASE places the box origin at 0,0,0 by default)
    lmp_trajectrories[k] = []
    for i,frame in enumerate(tmp_traj):
        frame.center() 
        lmp_trajectrories[k].append( frame )
    tmp_traj.close()
%memit
peak memory: 29625.96 MiB, increment: -4.30 MiB
In [110]:
%memit
peak memory: 29627.56 MiB, increment: 0.00 MiB
In [116]:
len(lmp_trajectrories['nvtEquilibration'])
Out[116]:
26
In [117]:
view(lmp_trajectrories['nvtEquilibration'], viewer='ase') # opens ASE GUI
In [118]:
#nv.show_asetraj(lmp_trajectrories['npt100ps'])
In [119]:
# Several problems with nglview:
#  1) does not display or infer bonding for ASE trajectory
#  2) even with displayed gui, not clear how to activate
In [120]:
lmp_trajectrories['npt1ns'][0]
Out[120]:
Atoms(symbols='C19H30452Au3024BrNO15205', pbc=True, cell=[63.027867814605706, 63.027867814605706, 126.05573562921141], masses=..., mol=...)

LAMMPS trajectories, stripped of solvent

In [121]:
%memit
# apparently, order of atoms changed in netcdf
lmp_naked_trajectrories = {}
for k in lmp_trajectrories:
    lmp_naked_trajectrories[k] = []
    for f in lmp_trajectrories[k]:
        #g = f.copy()
        #del g[water_indices]
        g = f[ (np.array(f.get_chemical_symbols()) != 'O') & (np.array(f.get_chemical_symbols()) != 'H')]
        lmp_naked_trajectrories[k].append(g)
peak memory: 29625.98 MiB, increment: 0.00 MiB
In [122]:
trajectoryView = nv.show_asetraj(lmp_naked_trajectrories['npt1ns'][::100])
In [123]:
#trajectoryView = nv.show_asetraj(lmp_naked_trajectrories['nvtEquilibration'][::100])
In [124]:
trajectoryView.remove_ball_and_stick()
In [125]:
trajectoryView.add_spacefill() # try a different representation sytle
In [126]:
trajectoryView
In [168]:
view(lmp_naked_trajectrories['npt1ns'], viewer='ase') # opens ASE GUI
In [169]:
view(lmp_naked_trajectrories['nvtEquilibration'], viewer='ngl') # opens ASE GUI
In [107]:
# Several problems with nglview:
#  1) does not display or infer bonding for ASE trajectory
#  2) even with displayed gui, not clear how to activate
#  3) somehow mixes up the atom type for netcdf

Make a movie via .png frames and ffmpeg

Preparation

In [127]:
# create a subdir from within the notebook
%mkdir png
mkdir: cannot create directory 'png': File exists
In [128]:
traj = lmp_naked_trajectrories['npt1ns']
In [129]:
nv.show_ase(traj[200])
In [130]:
# from ~ frame 3500 to the end ~ means
In [131]:
totalFramesAvailable = 5000
In [132]:
desiredVideoDuration = 30 # s
In [133]:
framesPerSecond = 30 # s^-1
In [134]:
neededFrames = desiredVideoDuration*framesPerSecond
In [135]:
neededFrames
Out[135]:
900
In [136]:
every_nth = np.ceil(totalFramesAvailable / neededFrames).astype(int)
In [137]:
every_nth
Out[137]:
6
In [140]:
png_prefix = absolute_prefix + os.sep + 'png' + os.sep + 'traj_1ns'

Orientation and bounding box settings by trial & error

In [148]:
f = traj[0].copy()
In [149]:
# find a desired orientation
#f.rotate('x', (-1,2,-1), rotate_cell=True)
#testframe.rotate(90, (1,0,0), rotate_cell=False)
#f.rotate('z', (1,1,-1), rotate_cell=False)
#f.rotate(-90,'x', rotate_cell=False)
# strange, somehow opposite behavior with "in-notebook" viewer and png renderer
In [150]:
f.rotate(-60,'x', rotate_cell=False)
In [151]:
f.rotate(30,'y', rotate_cell=False)
In [152]:
view(f,viewer='ngl')
In [389]:
nv.show_ase(f)
In [390]:
cell = f.get_cell()
In [391]:
cell
Out[391]:
array([[  63.02786781,    0.        ,    0.        ],
       [   0.        ,   63.02786781,    0.        ],
       [   0.        ,    0.        ,  126.05573563]])
In [158]:
bbox = [-20, 20, 80, 120 ]
In [159]:
# one trial
# the commented lines can be used to replicate the cell
#cell = f.get_cell()
#f = f.repeat((3, 1, 3))
#f.set_cell(cell)
# bbox = [-5, -5, cell[0,0] + 5, cell[1,1] + 5 ]
#f.center()
ase.io.write(png_prefix + '_test.png', f, show_unit_cell=False, 
                bbox=bbox)
# the bounding vox's 1st coordinate corresponds to the horizontal direction
# and ASE's x direction

Load image with

![title](png/traj_1ns_test.png?arg)

and change to some random string after the question mark in order to enforce reloading when image changed on disk ([https://github.com/jupyter/notebook/issues/1369]) title

Batch rendering

In [160]:
bbox = [-20, 20, 80, 120 ]
In [162]:
# make a movie
# https://wiki.fysik.dtu.dk/ase/development/making_movies.html
#for i,f in enumerate(frames_1ns_stripped[0::every_nth]):
for i,g in enumerate(traj[::every_nth]):
    #f.rotate('-y', 'z', rotate_cell=True)
    f = g.copy()
    f.rotate(-60,'x', rotate_cell=False)
    f.rotate(30,'y', rotate_cell=False)
    #f.center()

    #cell = f.get_cell()
    #f = f.repeat((3, 1, 3))
    #f.set_cell(cell)
    ase.io.write(png_prefix + '_{:05d}.png'.format(i), f, show_unit_cell=False, 
                bbox=bbox)
# bbox measures chosen to exceed cell once in x direction and twice in (repeated) y direction
# externally execute
#    ffmpeg -r 30 -f image2 -i "traj_1ns_%05d.png" -vcodec libx264 -crf 25 -pix_fmt yuv420p "traj_1ns.mp4
In [143]:
# operating some bash commands from within jupyter nb
%cd png
/work/ws/nemo/fr_jh1130-201708-0/jobs/lmplab/ctab/201807/1_CTAB_on_AU_21x12x2_netcdf/png
In [144]:
group_ws=%env GROUP_WS
In [ ]:
group_ws
In [ ]:
!find $group_ws -name ffmpeg
In [ ]:
%pwd
In [163]:
%%bash
source '/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17/env.sh'
cd /work/ws/nemo/fr_jh1130-201708-0/jobs/lmplab/ctab/201807/1_CTAB_on_AU_21x12x2_netcdf/png
ffmpeg -r 30 -f image2 -i "traj_1ns_%05d.png" -vcodec libx264 -crf 25 -pix_fmt yuv420p "traj_1ns.mp4"
# Setting up environment
# ----------------------
# Root directory:    /work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17
# Python executable: python3
# Python platform:   linux-x86_64-3.6
# PATH:              /work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17/bin
# LD_LIBRARY_PATH:   /work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17/lib
# CPATH:             /work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17/include
# LIBRARY_PATH:      /work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17/lib
# PYTHONPATH:        /work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17/lib/python3.6/site-packages
# Updating ~/.config/pip/pip.conf
ffmpeg version N-89659-g903be5e4f6 Copyright (c) 2000-2017 the FFmpeg developers
  built with gcc 5.2.0 (GCC)
  configuration: --prefix=/work/ws/nemo/fr_lp1029-IMTEK_SIMULATION-0/local_Nov17 --enable-libx264 --enable-shared --enable-gpl
  libavutil      56.  7.100 / 56.  7.100
  libavcodec     58.  9.100 / 58.  9.100
  libavformat    58.  3.100 / 58.  3.100
  libavdevice    58.  0.100 / 58.  0.100
  libavfilter     7.  8.100 /  7.  8.100
  libswscale      5.  0.101 /  5.  0.101
  libswresample   3.  0.101 /  3.  0.101
  libpostproc    55.  0.100 / 55.  0.100
Input #0, image2, from 'traj_1ns_%05d.png':
  Duration: 00:00:27.80, start: 0.000000, bitrate: N/A
    Stream #0:0: Video: png, rgba(pc), 2000x2000 [SAR 2834:2834 DAR 1:1], 30 fps, 30 tbr, 30 tbn, 30 tbc
Stream mapping:
  Stream #0:0 -> #0:0 (png (native) -> h264 (libx264))
Press [q] to stop, [?] for help
[libx264 @ 0x1de1ec0] using SAR=1/1
[libx264 @ 0x1de1ec0] using cpu capabilities: MMX2 SSE2Fast SSSE3 SSE4.2 AVX FMA3 BMI2 AVX2
[libx264 @ 0x1de1ec0] profile High, level 5.0
[libx264 @ 0x1de1ec0] 264 - core 155 r2893 b00bcaf - H.264/MPEG-4 AVC codec - Copyleft 2003-2017 - http://www.videolan.org/x264.html - options: cabac=1 ref=3 deblock=1:0:0 analyse=0x3:0x113 me=hex subme=7 psy=1 psy_rd=1.00:0.00 mixed_ref=1 me_range=16 chroma_me=1 trellis=1 8x8dct=1 cqm=0 deadzone=21,11 fast_pskip=1 chroma_qp_offset=-2 threads=30 lookahead_threads=5 sliced_threads=0 nr=0 decimate=1 interlaced=0 bluray_compat=0 constrained_intra=0 bframes=3 b_pyramid=2 b_adapt=1 b_bias=0 direct=1 weightb=1 open_gop=0 weightp=2 keyint=250 keyint_min=25 scenecut=40 intra_refresh=0 rc_lookahead=40 rc=crf mbtree=1 crf=25.0 qcomp=0.60 qpmin=0 qpmax=69 qpstep=4 ip_ratio=1.40 aq=1:1.00
Output #0, mp4, to 'traj_1ns.mp4':
  Metadata:
    encoder         : Lavf58.3.100
    Stream #0:0: Video: h264 (libx264) (avc1 / 0x31637661), yuv420p, 2000x2000 [SAR 1:1 DAR 1:1], q=-1--1, 30 fps, 15360 tbn, 30 tbc
    Metadata:
      encoder         : Lavc58.9.100 libx264
    Side data:
      cpb: bitrate max/min/avg: 0/0/0 buffer size: 0 vbv_delay: -1
frame=  834 fps= 34 q=-1.0 Lsize=   34593kB time=00:00:27.70 bitrate=10230.6kbits/s speed=1.12x    
video:34582kB audio:0kB subtitle:0kB other streams:0kB global headers:0kB muxing overhead: 0.031500%
[libx264 @ 0x1de1ec0] frame I:4     Avg QP:13.55  size:181543
[libx264 @ 0x1de1ec0] frame P:210   Avg QP:25.33  size: 61940
[libx264 @ 0x1de1ec0] frame B:620   Avg QP:29.63  size: 34964
[libx264 @ 0x1de1ec0] consecutive B-frames:  0.8%  0.0%  0.4% 98.8%
[libx264 @ 0x1de1ec0] mb I  I16..4: 66.9% 11.6% 21.5%
[libx264 @ 0x1de1ec0] mb P  I16..4:  0.5%  1.3%  0.5%  P16..4:  5.9%  9.6%  7.4%  0.0%  0.0%    skip:74.7%
[libx264 @ 0x1de1ec0] mb B  I16..4:  0.3%  0.4%  0.2%  B16..8:  8.2%  9.7%  4.9%  direct: 0.8%  skip:75.6%  L0:50.7% L1:40.8% BI: 8.5%
[libx264 @ 0x1de1ec0] 8x8 transform intra:39.9% inter:12.7%
[libx264 @ 0x1de1ec0] coded y,uvDC,uvAC intra: 15.4% 34.4% 27.3% inter: 5.6% 5.9% 5.1%
[libx264 @ 0x1de1ec0] i16 v,h,dc,p: 68% 16% 17%  0%
[libx264 @ 0x1de1ec0] i8 v,h,dc,ddl,ddr,vr,hd,vl,hu: 10%  7% 82%  0%  0%  0%  0%  0%  0%
[libx264 @ 0x1de1ec0] i4 v,h,dc,ddl,ddr,vr,hd,vl,hu: 19% 20% 32%  5%  5%  5%  6%  4%  5%
[libx264 @ 0x1de1ec0] i8c dc,h,v,p: 60% 21% 17%  1%
[libx264 @ 0x1de1ec0] Weighted P-Frames: Y:0.0% UV:0.0%
[libx264 @ 0x1de1ec0] ref P L0: 46.2%  4.0% 25.9% 23.9%
[libx264 @ 0x1de1ec0] ref B L0: 63.8% 26.1% 10.1%
[libx264 @ 0x1de1ec0] ref B L1: 81.1% 18.9%
[libx264 @ 0x1de1ec0] kb/s:10190.39

4. Radial distribution functions

4.0. Index selections

In [138]:
# helps in selecting atom nmbers based on element abbreviations
chem_sym = np.array(ase.data.chemical_symbols)
In [139]:
sulfur_index = [ a.number - 1 for a in surfactant if a.name == 'S' ][0]
In [140]:
tail_carbon_index = [ a.number - 1 for a in surfactant if a.name == 'C12' ][0]
In [141]:
sulfur_index
Out[141]:
144
In [142]:
tail_carbon_index
Out[142]:
182

4.1. Headgroup, Tailgroud - Gold RDF

In [143]:
# Quoted from https://wiki.fysik.dtu.dk/asap/Radial%20Distribution%20Functions
#
# Partial RDFs: Looking at specific elements or subsets of atoms
# It is often useful to look at partial RDFs, for example RDFs only taking some 
# elements into account, for example to get the distribution of atoms of element 
# B around an atom of element A. Do do this, call get_rdf() with the optional argument, 
# elements. It must be a tuple of two atomic numbers (a, b), the returned RDF then tells 
# how many b neighbors an a atom has.
#
# It is also possible to group atoms according to other criteria, for example to 
# calculate RDFs in different parts of space. In this case, the atoms must be divided 
# into groups when the RadialDistributionFunction object is created. Pass the extra 
# argument groups when creating the object, it must be an array containing a 
# non-negative integer for each atom specifying its group. When calling get_rdf() 
# use the argument group to specify for which group of atoms you want the RDF.

# IMPORTANT: The partial RDFs are normalized such that they sum up to the global RDF. 
# This means that integrating the first peak of a partial RDF obtained with 
# elements=(a,b) does not give you the number of B atoms in the first shell around 
# the A atom. Instead it gives this coordination number multiplied with the concentration 
# of A atoms.

traj = lmp_trajectrories['npt1ns']

# pick indices to look at during rdf computation
surfaceSoluteRdfIndices = [sulfur_index] +[tail_carbon_index] + surface_indices

nSegment = 1000

absoluteEnd = len(traj)-1
rdfSulfurGoldList     = []
rdfTailCarbonGoldList = []

# rMax is the rdf length (in Angstrom for LAMMPS output in real units)
rMax  = 20

# nBins: can be understood as the number of data points on the RDF
nBins = 1000
# actual distances
rdf_x = (np.arange(nBins) + 0.5) * rMax / nBins

# instead of computing an average rdf over the whole trajectory, 
# we split the trajectory into several timespans of nSegement timestep length
for nStart in range(0,absoluteEnd,nSegment):
    print(nStart) # some progress report
    surfaceSoluteRdf = None
    for frame in traj[nStart:(nStart+nSegment)]:
        # the asap rdf functionality is not that convenient, but explicitely choosing
        # only the atoms we are interested in a priori, we can get exactly the rdf we want by 
        # making use of the "elements" option
        if surfaceSoluteRdf is None:
            surfaceSoluteRdf = RadialDistributionFunction(frame[surfaceSoluteRdfIndices],
                                     rMax = rMax, nBins = nBins)
        else:
            surfaceSoluteRdf.atoms = frame[surfaceSoluteRdfIndices]  # Fool RDFobj to use the new atoms
        surfaceSoluteRdf.update()           # Collect data
    
    # np.where facilitates the selection of according atom numbers by specifying the chemical symbol
    rdfSulfurGold = surfaceSoluteRdf.get_rdf(elements=(
        np.where(chem_sym == 'S')[0][0],
        np.where(chem_sym == 'Au')[0][0] ))
    rdfTailCarbonGold = surfaceSoluteRdf.get_rdf(elements=(
        np.where(chem_sym == 'C')[0][0],
        np.where(chem_sym == 'Au')[0][0] ))
    
    rdfSulfurGoldList.append(rdfSulfurGold)
    rdfTailCarbonGoldList.append(rdfTailCarbonGold)
0
1000
2000
3000
4000
In [176]:
# code snippet for neat plotting of all time-segemtn rdfs
cols = 2
rows = np.ceil(len(rdfSulfurGoldList)/cols).astype(int)
pos = subplotPosition(rows,cols)
fig = plt.figure(figsize=(rows*5,cols*8))

for i, (rdfSulfurGold, rdfTailCarbonGold) in enumerate(zip(rdfSulfurGoldList,rdfTailCarbonGoldList)):
    p = next(pos)
    _, ax = addSubplot(rdf_x,rdfSulfurGold, 
                       legend = "S - Au RDF", 
                       xlabel = r'$\frac{r}{\AA}$', 
                       ylabel='arbitrary density', 
                       title = "{} ps - {} ps".format(i*nSegment,(i+1)*nSegment),
                       fig = fig, pos = p)
    _, _ = addSubplot(rdf_x, rdfTailCarbonGold, 
                      legend="S - tail C RDF",
                      ax = ax, pos = p)
    
fig.tight_layout()

4.2. head group sulfur, tail group carbon - water RDF

In [177]:
# Element tuples
element_tuples = [
    ( np.where(chem_sym == 'S')[0][0], np.where(chem_sym == 'O')[0][0] ),
    ( np.where(chem_sym == 'C')[0][0], np.where(chem_sym == 'O')[0][0] ) ]
In [178]:
element_tuples # in atomic numbers
Out[178]:
[(16, 8), (6, 8)]
In [179]:
surfactantSolventIndicesOfInterest = water_indices + [ sulfur_index ] + [ tail_carbon_index ]
In [180]:
surfactantSolventRDFs, surfactantSolventRDFx, surfactantSolventRDFobj = piecewiseRDF(
    lmp_trajectrories['npt1ns'], surfactantSolventIndicesOfInterest, element_tuples)
0
1000
2000
3000
4000
In [181]:
plotPiecewiceRdf(surfactantSolventRDFx, surfactantSolventRDFs, 
                 legend= [ "head group sulfur - water RDF", "tail group carbon - water RDF"]);

4.3. Sodium counterion RDF

In [182]:
# Element tuples
counterionRdf_element_tuples = [
    ( np.where(chem_sym == 'Na')[0][0], np.where(chem_sym == 'S')[0][0] ),
    ( np.where(chem_sym == 'Na')[0][0], np.where(chem_sym == 'C')[0][0] ),
    ( np.where(chem_sym == 'Na')[0][0], np.where(chem_sym == 'O')[0][0] ),
    ( np.where(chem_sym == 'Na')[0][0], np.where(chem_sym == 'Au')[0][0] ) ]
In [183]:
counterionRdf_element_tuples
Out[183]:
[(11, 16), (11, 6), (11, 8), (11, 79)]
In [184]:
counterionRdfLabels = [ "Na+ counterion - SDS head sulfur RDF", 
                        "Na+ counterion - SDS tail carbon RDF",
                        "Na+ counterion - water oxygen RDF",
                        "Na+ counterion - surface gold RDF"]
In [185]:
counterionRdfIndicesOfInterest = water_indices + surface_indices + \
    ion_indices + [ sulfur_index ] + [ tail_carbon_index ]
In [186]:
counterionRDFs, counterionRDFx, counterionRDFobj = piecewiseRDF(
    lmp_trajectrories['npt1ns'], counterionRdfIndicesOfInterest, counterionRdf_element_tuples)
0
1000
2000
3000
4000
In [187]:
plotPiecewiceRdf(counterionRDFx, counterionRDFs, legend= counterionRdfLabels);
In [189]:
plotPiecewiceRdf(counterionRDFx, counterionRDFs[0:2], 
                 legend= (counterionRdfLabels[0:2]));
In [190]:
plotPiecewiceRdf(counterionRDFx, counterionRDFs[2:], 
                 legend= (counterionRdfLabels[2:]));

5. Distance analysis

5.1. Headgroup - gold distance

Head approaches surfaces, apparently "stepwise"

In [191]:
traj = lmp_trajectrories['npt1ns']
In [192]:
traj[0][surfactant_indices].get_atomic_numbers()
Out[192]:
array([16,  8,  8, ...,  1,  1,  1])
In [193]:
averageDistanceS2Au, averageDistanceS2AuTimes = piecewiseAveragedDistance(traj,
                                reference_index=sulfur_index,
                                atom_indices=surface_indices,
                                nSegment=50)
In [194]:
len(averageDistanceS2AuTimes)
Out[194]:
100
In [195]:
averageDistanceS2Au.shape
Out[195]:
(3, 100)
In [196]:
distanceLabels = ['x', 'y', 'z']
In [197]:
distanceFig = plt.figure()
for i in range(0,3):
    plt.plot( averageDistanceS2AuTimes*0.2, averageDistanceS2Au[i,:], label= distanceLabels[i] )
plt.legend()
plt.title("Distance of head group sulfur from gold atoms")
plt.xlabel("time t / ps")
plt.ylabel(r'Distance $\frac{r}{\AA}$')
Out[197]:
Text(0,0.5,'Distance $\\frac{r}{\\AA}$')
In [198]:
traj = lmp_trajectrories['npt1ns']
In [199]:
traj[0][surfactant_indices].get_atomic_numbers()
Out[199]:
array([16,  8,  8, ...,  1,  1,  1])
In [201]:
# slight decrease in potential in comparison with approach towards surface
# confirms the anticipated energetically favored adsorption state
nptProduction_1ns_thermo_pd[["PotEng","E_pair"]].rolling(window=5000,center=True).mean().plot()
Out[201]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fd03057f6d8>
In [202]:
nptProduction_1ns_thermo_pd
Out[202]:
TotEng KinEng PotEng Temp Press Enthalpy E_bond E_angle E_dihed E_impro E_pair E_vdwl E_coul E_long E_tail Volume E_intramolecular
Step
0 -6263.0571 453.10079 -6716.1579 301.79866 798.122940 -6171.1625 2.914661 35.266749 7.031867 0 -6761.3711 -4897.9271 7191.5053 -9054.9493 0 7894.8534 45.213277
10 -6265.7761 453.21357 -6718.9897 301.87378 -119.442190 -6279.6182 4.618620 31.398426 7.081973 0 -6762.0887 -4903.4927 7199.2171 -9057.8131 0 7946.3381 43.099019
20 -6263.2673 460.64210 -6723.9094 306.82173 289.621950 -6229.7871 2.526166 35.519864 8.142813 0 -6770.0983 -4897.2971 7183.2694 -9056.0706 0 7926.4993 46.188843
30 -6267.9400 458.45626 -6726.3962 305.36580 -1420.633000 -6433.0946 3.700661 32.282913 8.018788 0 -6770.3986 -4912.9657 7200.8081 -9058.2410 0 7971.3700 44.002362
40 -6267.0365 454.66485 -6721.7013 302.84044 1557.912000 -6086.8862 2.867427 30.681025 8.562313 0 -6763.8121 -4884.7728 7175.9668 -9055.0061 0 7928.9562 42.110766
50 -6267.5204 446.52921 -6714.0497 297.42150 -197.479070 -6290.4505 2.329377 36.206542 8.229080 0 -6760.8147 -4888.3573 7184.5147 -9056.9720 0 7961.7322 46.765000
60 -6274.3052 436.75049 -6711.0557 290.90815 -1806.958100 -6485.3207 2.942238 32.878879 8.945283 0 -6755.8221 -4879.5523 7179.8467 -9056.1165 0 8007.3771 44.766400
70 -6274.5247 432.79404 -6707.3187 288.27286 2262.071700 -6012.4026 3.002364 31.435201 7.462704 0 -6749.2190 -4839.6187 7147.0856 -9056.6858 0 7945.5032 41.900269
80 -6271.3038 432.60558 -6703.9094 288.14733 -767.276420 -6360.7215 3.429762 35.319369 8.555790 0 -6751.2143 -4859.1891 7165.5075 -9057.5327 0 7990.8945 47.304921
90 -6277.3733 436.87190 -6714.2452 290.98902 -1400.424600 -6440.4945 4.859236 36.321427 8.025444 0 -6763.4513 -4872.0685 7165.4209 -9056.8037 0 7986.8356 49.206106
100 -6271.7345 461.29208 -6733.0265 307.25467 2055.301600 -6034.5404 2.256737 33.865517 7.711067 0 -6776.8599 -4884.5035 7164.3725 -9056.7288 0 7913.2039 43.833321
110 -6273.5742 463.61085 -6737.1851 308.79914 -323.591330 -6311.0655 2.352985 28.726380 8.388930 0 -6776.6534 -4906.8128 7185.2038 -9055.0443 0 7944.3293 39.468295
120 -6276.0452 461.65982 -6737.7050 307.49961 -441.847120 -6327.1367 5.778075 33.019169 8.394795 0 -6784.8970 -4888.5805 7159.4262 -9055.7427 0 7928.6795 47.192039
130 -6280.6749 459.68075 -6740.3557 306.18140 -637.112930 -6354.4156 6.401045 30.353133 8.444789 0 -6785.5546 -4883.2984 7152.9973 -9055.2535 0 7936.2409 45.198967
140 -6283.6515 460.93064 -6744.5821 307.01392 1156.133700 -6150.3348 4.649919 35.009942 8.537100 0 -6792.7791 -4872.3258 7137.1999 -9057.6533 0 7906.7986 48.196961
150 -6287.4721 454.13376 -6741.6059 302.48669 -756.337410 -6374.9679 2.203549 39.599300 8.247860 0 -6791.6566 -4884.3972 7149.9175 -9057.1769 0 7932.2335 50.050709
160 -6296.9662 452.95837 -6749.9246 301.70380 -801.308860 -6389.5545 2.023112 36.966138 7.017698 0 -6795.9315 -4897.3117 7158.7319 -9057.3517 0 7922.8241 46.006948
170 -6298.2876 430.20906 -6728.4966 286.55107 1931.707400 -6076.7489 3.677957 34.455475 7.294375 0 -6773.9244 -4879.2148 7161.7515 -9056.4611 0 7863.7982 45.427807
180 -6302.9410 437.23602 -6740.1770 291.23155 -1420.521800 -6466.9540 4.149034 35.587254 7.881505 0 -6787.7948 -4905.1402 7176.5710 -9059.2255 0 7916.8895 47.617793
190 -6307.0221 431.59297 -6738.6150 287.47286 915.942100 -6202.1984 3.741286 33.683479 10.019462 0 -6786.0593 -4904.5362 7172.4806 -9054.0037 0 7847.2127 47.444227
200 -6306.9505 430.02666 -6736.9772 286.42958 -1258.643600 -6452.1378 4.410197 33.476811 9.002909 0 -6783.8671 -4905.8964 7180.3927 -9058.3633 0 7909.5128 46.889917
210 -6309.5913 420.09868 -6729.6900 279.81681 1169.287600 -6175.6618 2.922126 31.358654 8.352171 0 -6772.3229 -4883.1658 7165.2451 -9054.4022 0 7853.7877 42.632951
220 -6303.3100 421.50866 -6724.8187 280.75596 72.554129 -6294.9682 2.952090 30.625387 9.300480 0 -6767.6966 -4869.3136 7158.8155 -9057.1986 0 7883.6208 42.877957
230 -6303.3970 422.26880 -6725.6658 281.26228 -1197.180300 -6441.3359 4.219537 37.778284 8.775535 0 -6776.4392 -4860.3210 7140.7967 -9056.9149 0 7900.4366 50.773356
240 -6293.9270 423.41199 -6717.3390 282.02372 1159.177500 -6161.0547 5.194174 38.834911 8.837476 0 -6770.2055 -4854.4758 7140.3737 -9056.1034 0 7859.7484 52.866561
250 -6287.0163 441.90500 -6728.9213 294.34143 -622.037030 -6358.6038 2.230116 41.752217 8.881751 0 -6781.7853 -4866.1292 7142.1529 -9057.8091 0 7891.2436 52.864084
260 -6278.9437 449.59318 -6728.5368 299.46234 4.150381 -6278.4675 3.769346 33.358466 8.309849 0 -6773.9745 -4871.3214 7154.7832 -9057.4363 0 7867.2950 45.437661
270 -6270.4279 438.30711 -6708.7350 291.94498 235.898650 -6243.3468 2.746512 38.512063 8.062146 0 -6758.0557 -4883.4178 7181.7441 -9056.3820 0 7871.6379 49.320721
280 -6260.3922 430.99919 -6691.3914 287.07736 265.349570 -6229.8798 3.209129 38.922099 8.821513 0 -6742.3441 -4873.4247 7185.7386 -9054.6580 0 7884.6529 50.952741
290 -6251.5076 442.33183 -6693.8394 294.62574 40.789692 -6246.8122 2.076750 36.608942 9.740352 0 -6742.2654 -4866.2202 7178.6422 -9054.6874 0 7892.9944 48.426044
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
499710 -6253.0134 437.40051 -6690.4139 291.34111 741.857250 -6165.8642 4.495545 38.366228 7.745658 0 -6741.0213 -4915.8919 7229.4598 -9054.5892 0 8055.0307 50.607431
499720 -6254.8108 459.86198 -6714.6728 306.30211 181.504630 -6233.5405 3.749332 43.081015 7.654077 0 -6769.1572 -4919.5434 7206.1765 -9055.7903 0 8035.4547 54.484424
499730 -6259.0044 462.36861 -6721.3730 307.97171 -1675.142400 -6456.5837 5.617976 34.245373 8.365765 0 -6769.6021 -4924.1408 7209.1965 -9054.6578 0 8087.4935 48.229114
499740 -6262.5737 450.90586 -6713.4796 300.33667 2559.336000 -5963.9443 5.006742 36.079960 7.870453 0 -6762.4367 -4892.6536 7185.5991 -9055.3822 0 8000.7258 48.957155
499750 -6263.8812 440.69591 -6704.5771 293.53609 -908.432570 -6371.0167 6.431512 38.878450 8.785731 0 -6758.6728 -4899.6812 7198.5829 -9057.5745 0 8086.5742 54.095693
499760 -6270.7581 426.86339 -6697.6215 284.32261 -494.337360 -6328.8817 3.836806 40.432847 7.816843 0 -6749.7080 -4891.4959 7197.6188 -9055.8309 0 8062.1940 52.086495
499770 -6271.5294 447.08539 -6718.6148 297.79196 546.610540 -6207.3319 3.413413 34.671077 8.060081 0 -6764.7594 -4907.4593 7199.1281 -9056.4282 0 8053.1259 46.144571
499780 -6272.5232 445.89931 -6718.4225 297.00194 850.320420 -6172.9345 6.490952 35.483908 8.762540 0 -6769.1599 -4902.5074 7189.6665 -9056.3189 0 8030.6650 50.737399
499790 -6274.3266 452.19365 -6726.5203 301.19444 -2031.606600 -6513.6513 5.695687 30.944562 7.674509 0 -6770.8351 -4910.5201 7194.5744 -9054.8894 0 8077.4070 44.314758
499800 -6279.9190 437.87174 -6717.7908 291.65499 14.469685 -6278.2220 2.733233 36.153346 7.932595 0 -6764.6099 -4896.6472 7189.6687 -9057.6315 0 8041.9307 46.819174
499810 -6273.7231 426.25198 -6699.9751 283.91536 1710.929800 -6073.9294 4.421265 44.999394 7.093955 0 -6756.4897 -4886.3751 7182.3848 -9052.4994 0 8007.0694 56.514614
499820 -6275.1466 434.14931 -6709.2959 289.17557 -683.052470 -6355.2914 4.538397 36.801849 7.709994 0 -6758.3461 -4897.9921 7195.4448 -9055.7988 0 8045.3626 49.050240
499830 -6272.3457 446.22269 -6718.5684 297.21734 340.794890 -6232.5669 3.601451 41.862426 7.867179 0 -6771.8995 -4894.5539 7177.8946 -9055.2402 0 8003.5582 53.331056
499840 -6272.6501 432.20283 -6704.8529 287.87907 -856.983810 -6373.1155 4.670558 42.418385 7.270109 0 -6759.2120 -4907.8847 7205.7972 -9057.1245 0 8038.3725 54.359052
499850 -6266.1940 444.33626 -6710.5303 295.96083 2374.045800 -5990.7949 6.065900 38.592813 8.181420 0 -6763.3704 -4880.0353 7171.6915 -9055.0267 0 7954.2200 52.840133
499860 -6263.2341 441.11070 -6704.3448 293.81237 -2625.696400 -6570.8756 2.850867 36.681695 7.316487 0 -6751.1939 -4920.3032 7222.7747 -9053.6654 0 8033.8637 46.849049
499870 -6264.5047 459.49249 -6723.9971 306.05600 913.565600 -6158.6473 4.996637 40.930141 9.431102 0 -6779.3550 -4907.7904 7183.1123 -9054.6769 0 7945.2132 55.357880
499880 -6255.2834 449.11794 -6704.4013 299.14579 384.137790 -6210.5924 2.533438 38.575715 8.237135 0 -6753.7476 -4899.3004 7199.3852 -9053.8324 0 7977.3182 49.346288
499890 -6254.2658 456.27998 -6710.5458 303.91623 -697.700950 -6335.3698 8.011759 38.542578 10.909320 0 -6768.0094 -4918.7211 7204.8636 -9054.1520 0 7970.7134 57.463657
499900 -6252.9597 447.24667 -6700.2063 297.89938 17.087292 -6250.9691 4.025048 37.090560 9.948318 0 -6751.2703 -4907.9099 7210.1227 -9053.4830 0 7987.9642 51.063925
499910 -6247.4622 443.81482 -6691.2771 295.61352 1603.568000 -6061.3013 4.503919 39.164993 9.094639 0 -6744.0406 -4897.0974 7208.1188 -9055.0620 0 7960.2257 52.763551
499920 -6245.7503 476.56218 -6722.3125 317.42568 -1342.857700 -6402.5542 4.722325 39.548960 7.774728 0 -6774.3585 -4914.6904 7195.8020 -9055.4701 0 8006.6505 52.046013
499930 -6250.0141 467.92444 -6717.9386 311.67231 -531.017620 -6312.0684 3.541105 40.314543 7.888704 0 -6769.6829 -4905.4988 7191.1729 -9055.3571 0 8012.8458 51.744352
499940 -6246.0720 477.54644 -6723.6184 318.08128 1897.185400 -6025.1890 7.547013 33.760970 8.207800 0 -6773.1342 -4898.4884 7180.8605 -9055.5063 0 7983.1905 49.515783
499950 -6250.9215 470.11093 -6721.0324 313.12867 -1863.994100 -6470.1134 4.184154 35.056179 7.937175 0 -6768.2099 -4922.4589 7208.8954 -9054.6465 0 8063.1365 47.177508
499960 -6257.8931 451.16495 -6709.0581 300.50925 237.301890 -6230.1499 4.812108 42.722978 8.651148 0 -6765.2443 -4927.4951 7218.7608 -9056.5101 0 8016.4116 56.186234
499970 -6260.3941 446.92707 -6707.3212 297.68650 564.970480 -6194.2917 3.059662 40.446686 11.477817 0 -6762.3054 -4941.8806 7235.3359 -9055.7606 0 8022.6104 54.984166
499980 -6267.3801 442.31905 -6709.6991 294.61723 756.080360 -6178.9228 5.148024 34.304925 10.712310 0 -6759.8644 -4938.8116 7232.7306 -9053.7833 0 8022.1318 50.165259
499990 -6269.9960 441.97783 -6711.9739 294.38994 -362.199670 -6312.4936 4.984028 33.503398 9.790720 0 -6760.2520 -4933.1726 7227.7669 -9054.8463 0 8045.2582 48.278147
500000 -6279.8060 448.26221 -6728.0682 298.57581 -745.332870 -6367.3912 4.698220 32.922787 8.444214 0 -6774.1335 -4919.4227 7200.2035 -9054.9143 0 8057.5780 46.065221

50001 rows × 17 columns

In [203]:
nptProduction_1ns_thermo_pd[["E_intramolecular"]].rolling(window=5000,center=True).mean().plot()
Out[203]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fd0510d0da0>
In [204]:
ax = nptProduction_1ns_thermo_pd[["E_intramolecular"]].plot()
nptProduction_1ns_thermo_pd[["E_intramolecular"]].rolling(window=5000,center=True).mean().plot(ax=ax)
Out[204]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fd06811bbe0>
In [205]:
# running average, non-bonded energy
ax = nptProduction_1ns_thermo_pd[["E_pair"]].plot()
nptProduction_1ns_thermo_pd[["E_pair"]].rolling(window=10000,center=True).mean().plot(ax=ax)
Out[205]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fd02349f208>
In [206]:
# running average, total energy
ax = nptProduction_1ns_thermo_pd[["TotEng"]].plot()
nptProduction_1ns_thermo_pd[["TotEng"]].rolling(window=10000,center=True).mean().plot(ax=ax)
Out[206]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fd023e3d780>
In [207]:
# running average, total energy
ax = nptProduction_1ns_thermo_pd[["PotEng"]].plot()
nptProduction_1ns_thermo_pd[["PotEng"]].rolling(window=10000,center=True).mean().plot(ax=ax)
Out[207]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fd01d0211d0>

5.2. Tailgroup - gold distance

z - direction: does not change much

In [210]:
averageDistanceTailC2Au, averageDistanceTailC2AuTimes = piecewiseAveragedDistance(traj,
                                reference_index=tail_carbon_index,
                                atom_indices=surface_indices,
                                nSegment=50)
In [211]:
for i in range(0,3):
    plt.plot( averageDistanceTailC2AuTimes*0.2, averageDistanceTailC2Au[i,:], label= distanceLabels[i] )
plt.legend()
plt.title("Distance of tail group from gold atoms")
plt.xlabel("time t / ps")
plt.ylabel(r'Distance $\frac{r}{\AA}$')
Out[211]:
Text(0,0.5,'Distance $\\frac{r}{\\AA}$')

5.3. Surface COM - Surfactant COM distance

In [212]:
averageDistanceComCom, averageDistanceComComTimes = piecewiseAveragedComComDistance(traj,
                                surfactant_indices,surface_indices,
                                nSegment=50)
In [213]:
for i in range(0,3):
    plt.plot( averageDistanceComComTimes*0.2, averageDistanceComCom[i,:], label= distanceLabels[i] )
plt.legend()
plt.title("Distance of SDS center of mass from gold layer center of mass")
plt.xlabel("time t / ps")
plt.ylabel(r'Distance $\frac{r}{\AA}$')
Out[213]:
Text(0,0.5,'Distance $\\frac{r}{\\AA}$')
In [215]:
# for comparison: Headgroup - surface distance: 
distanceFig
Out[215]:

SDS chain length

In [216]:
averageChainLength, averageChainLengthTimes = piecewiseAveragedDistance(traj,
                                reference_index=tail_carbon_index,
                                atom_indices=[sulfur_index],
                                nSegment=50)
In [217]:
np.linalg.norm(averageChainLength,axis=0)
Out[217]:
array([ 13.19284389,  13.30118302,  11.45347917, ...,  15.15614958,
        15.05356845,  15.03744343])
In [218]:
np.linalg.norm(averageChainLength,axis=0).shape
Out[218]:
(100,)
In [219]:
plt.plot( averageChainLengthTimes*0.2, 
         np.linalg.norm(averageChainLength,axis=0), label="SDS Chain length")
plt.legend()
plt.title("Distance of SDS head group sulfur to tail group carbon")
plt.xlabel("time t / ps")
plt.ylabel(r'Distance $\frac{r}{\AA}$')
Out[219]:
Text(0,0.5,'Distance $\\frac{r}{\\AA}$')

MSD and diffusivities

In [220]:
dt
Out[220]:
2e-15
In [221]:
T = 1e-9 # 1 ns
In [222]:
Nf = len(traj) - 1 # number of stored frames, corresponds to 1ns
In [223]:
Nf
Out[223]:
5000
In [224]:
Ns = T/dt # number of steps
In [225]:
Ns
Out[225]:
500000.0
In [226]:
StepsPerFrame = Ns / Nf
In [227]:
StepsPerFrame
Out[227]:
100.0
In [228]:
TimePerFrame = StepsPerFrame*dt
In [229]:
TimePerFrame
Out[229]:
2e-13
In [235]:
# displacements over fixed time spans, i.e. from each frame to the 50th or 5th following frame
In [230]:
displacement10ps = comDisplacement(traj, surfactant_indices, dt=50)
In [231]:
displacement1ps = comDisplacement(traj, surfactant_indices, dt=5)
In [232]:
# instead of averaging over many trajectories, we average over N time-wise close frames-frame displacements
# crude (all displacements within each set very correlated)
# default N = 500 ~ 1 ps
evaluateDisplacement(displacement10ps, dt=50);
# however, the results give a good idea about anisotropic and decreasing mobility 
# during approach towards surface
In [233]:
evaluateDisplacement(displacement1ps, dt=5);
In [234]:
evaluateDisplacement(displacement1ps, dt=5, window=1000);